Mal_a
Mal_a

Reputation: 3760

R: The difference between two points from the plot for log x axis

I need to calculate the time difference of the each sample (the samples can be differentiated using ID column) between two temperature points for log x axis. I get the calculated results, however for linear axis, not log. How can I achieve calculation for log scale using my code:

  1. Sample data

    dput(data)
    structure(list(id = c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,2L, 2L),
    Zeit = c(0L, 180L, 360L, 420L, 600L, 604L, 0L, 180L,360L, 480L, 600L,
    605L), Temp = c(963L, 824L, 666L, 658L, 641L,549L, 957L, 823L, 661L,
    660L, 642L, 562L)), .Names = c("id","Zeit", "Temp"), row.names = c(NA,
    12L), class = "data.frame")
    
  2. Code:

    Zt <- vapply(unique(data$id), function(ID){
        with(data[data$id == ID,], approx(x = Temp, y = Zeit, xout = 600))$y
    }, double(1))
    data.frame(id = unique(data$id), time = Zt)
    

There is an option in approx for method. However there are only two methods specified:

and as I mentioned before, the log is what I am looking for

Thanks for help!

[UPDATE]

Why log makes difference in my case.

Just for explanation I am going to use small set of data where we can see how log makes the difference:

here is the data:

data <-structure(list(id = c(1L, 1L, 1L), Zeit = c(31L, 701L, 902L), 
    Temp = c(930L, 549L, 481L)), .Names = c("id", "Zeit", "Temp"
), row.names = c(NA, 3L), class = "data.frame")

Here is plot for it (1st is with log axis, second is normal):

enter image description here enter image description here

if i wanna see for example the time which i need to get to 700°C, in log scale it would be equal to around 200 sec, for normal one it is equal to around 325 sec.

Upvotes: 1

Views: 333

Answers (2)

joemienko
joemienko

Reputation: 2280

A couple of preliminary items:

  1. Based on the graphs you are displaying below, your x and y values appear to be reversed in your original code.

  2. The function I have here does not have all of the functionality of the approx() function in stats, but I think it will meet your needs.

To start, consider that the approx() function uses a variant of the following formula to produce a result:

Equation 1:

where

To return a value of http://latex.codecogs.com/gif.download?%24y_%7Bout%7D%24 which corresponds to a log-scaled x-axis, we simply log the relevant portions of the formula as follows:

Equation 2:

Below I implement both formulas in a new function called approx_log()

approx_log <- function(x, y, xout){
  
  dat <- data.frame(y=y, x=x)
  dat <- dat[with(dat, order(x, y)), ]
  
  y_in <- dat$y
  x_in <- dat$x
  
  # find the start of our interval
  int_start <- which(x_in == max(x_in[x_in <= xout]))

  # assign the int_start value to x_0 and the 
  # value from the next highest index to x_1
  x_0 <- x_in[int_start]
  x_1 <- x_in[int_start + 1]
  
  # repeat for corresponding y-values 
  y_0 <- y_in[int_start]
  y_1 <- y_in[int_start + 1]
  

  y_out_lin <- y_0 + ((xout-x_0)/(x_1-x_0))*(y_1-y_0)
    
  y_out_log <- y_0 + ((log(xout)-log(x_0))/(log(x_1)-log(x_0)))*(y_1-y_0)
    
  # return values 
  list(x = xout, y_lin = y_out_lin, y_log = y_out_log)
}

As can be seen, this function returns a list of http://latex.codecogs.com/gif.download?x_%7Bout%7D_ and the log and linear interpolation values of http://latex.codecogs.com/gif.download?%24y_%7Bout%7D%24. The code below tests the function based on your visual interpolation in the posting.

data <-structure(list(id = c(1L, 1L, 1L), Zeit = c(31L, 701L, 902L), 
    Temp = c(930L, 549L, 481L)), .Names = c("id", "Zeit", "Temp"
), row.names = c(NA, 3L), class = "data.frame")

approx_log(x = data$Zeit, y = data$Temp, xout = 200)

## $x
## [1] 200
## 
## $y_lin
## [1] 833.897
## 
## $y_log
## [1] 702.2286

As you suggest, the log-scaled version of time at 200 seconds corresponds almost exactly to 700 degrees.

approx_log(x = data$Zeit, y = data$Temp, xout = 325)

## $x
## [1] 325
## 
## $y_lin
## [1] 762.8149
## 
## $y_log
## [1] 642.9125

The linear version of time at 325 seconds is a little higher (about 763 degrees) but reasonable based on your original plot. As a sanity check, we can see that the linear value matches the approx() function exactly.

approx(x = data$Zeit, y = data$Temp, xout = 325)

## $x
## [1] 325
## 
## $y
## [1] 762.8149

We can also run this through vapply() per your original request.

data <- structure(list(id = c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,2L, 2L),
               Zeit = c(0L, 180L, 360L, 420L, 600L, 604L, 0L, 180L,360L, 480L, 600L,
                        605L), Temp = c(963L, 824L, 666L, 658L, 641L,549L, 957L, 823L, 661L,
                                        660L, 642L, 562L)), .Names = c("id","Zeit", "Temp"), row.names = c(NA,
                                                                                                           12L), class = "data.frame")

Zt <- vapply(unique(data$id), function(ID){
  with(data[data$id == ID,], approx_log(y = Temp, x = Zeit, xout = 325))$y_log
}, double(1))
data.frame(id = unique(data$id), time = Zt)

##   id     time
## 1  1 689.3140
## 2  2 684.9043

You could seperately extract the y_lin value for the sake of comparison.

Zt <- vapply(unique(data$id), function(ID){
  with(data[data$id == ID,], approx_log(y = Temp, x = Zeit, xout = 325))$y_lin
}, double(1))
data.frame(id = unique(data$id), time = Zt)

##   id     time
## 1  1 696.7222
## 2  2 692.5000

EDIT The original question sought to solve for Zeit, given Temp (i.e. solve for x, given y). The above code solves for log-interpolated values of y, for a given x. The inverse of this is achieved by re-arranging Equation 2 to solve for x_out. The function for this, as provided in a comment by @joemienko is

x_out_log <- exp((yout*log(x_0)-y_1*log(x_0)-yout*log(x_1)+y_0*log(x_1))/(y_0-y_1))

Upvotes: 3

dww
dww

Reputation: 31452

The simplest way to do this is to use approx(). But you need to do a couple of things to get the behaviour you want:

  1. log transform the input values before passing them to approx, and then
  2. inverse-log (i.e. exp) the answer to convert back to the correct units.

Thus, in your example, to solve for Zeit at Temp=700, by interpolating the log-transformed values of Zeit, this would be

exp( approx(x = data$Temp, y = log(data$Zeit), xout = 700))$y )
## 203.6818

Upvotes: 1

Related Questions