[R] Error running caret's gbm train function with new version of caret
I am running caret for model exploration. I developed my code a number of
months ago and I've been running it with no issues. Recently, I updated my
version of caret however, and now I am getting a new error. I'm wondering
if this is due to the new release.
The error I am getting is when I am running GBM.
print(paste("calculating GBM for", i))
#gbm runs over and over again
set.seed(1)
trainModelGBM <- train(trainClass3, trainAsym, "gbm", metric="RMSE",
tuneLength = 5, trControl = con)
The error I am getting is at the end of the run once all the iterations
have been processed:
Error in { :
task 1 failed - "arguments imply differing number of rows: 5, 121"
trainClass3 and trainAsym have 311 values in them. I'm using 5 variables in
my matrix. I'm not sure where the 117 is coming from.
I found solutions online that suggested that updated the version of glmnet,
Matrix and doing something with cv.folds would work. None of these
solutions have worked for me.
Here is my R session info.
R version 2.15.1 (2012-06-22)
Platform: x86_64-unknown-linux-gnu (64-bit)
caret version 5.15-61
Thank you,
Katrina
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Installing rgdal in R: correct -configure flags for GDAL install on Linux Redhat
Hi, I'm installing rgdal but I keep having failures because I have not been
able to find a good source of information for the correct configuration
settings when installing GDAL.
My error from the R install.packages("rgdal") is below.
Can someone point me to a good source to tell me how to set after
./configure when installing GDAL?
I'd like to be able to work with raster images, geotiffs, netCDF files, and
other raster-based image processing in R.
Right now I'm running ./configure for GDAL as follows:
./configure --prefix=$HOME/local/gdal/gdal-1.8.0 --with
jasper=$HOME/local/jasper/jasper-1.900.1.uuid/src/libjasper
Any insight would be greatly appreciated! Thank you.
Error output from R:
* installing *source* package 'rgdal' ...
gdal-config: gdal-config
checking for gcc... gcc -std=gnu99
checking for C compiler default output file name... a.out
checking whether the C compiler works... yes
checking whether we are cross compiling... no
checking for suffix of executables...
checking for suffix of object files... o
checking whether we are using the GNU C compiler... yes
checking whether gcc -std=gnu99 accepts -g... yes
checking for gcc -std=gnu99 option to accept ANSI C... none needed
checking how to run the C preprocessor... gcc -std=gnu99 -E
checking for egrep... grep -E
checking for ANSI C header files... yes
checking for sys/types.h... yes
checking for sys/stat.h... yes
checking for stdlib.h... yes
checking for string.h... yes
checking for memory.h... yes
checking for strings.h... yes
checking for inttypes.h... yes
checking for stdint.h... yes
checking for unistd.h... yes
checking proj_api.h usability... no
checking proj_api.h presence... no
checking for proj_api.h... no
Error: proj_api.h not found.
If the PROJ.4 library is installed in a non-standard location,
use --configure-args='--with-proj-
include=/opt/local/include'
for example, replacing /opt/local/* with appropriate values
for your installation. If PROJ.4 is not installed, install it.
ERROR: configuration failed for package 'rgdal'
* removing
'/import/home/u1/uaf/kbennett/R/x86_64-unknown-linux-gnu-library/2.11/rgdal'
--
Katrina E. Bennett
PhD Student
University of Alaska Fairbanks
International Arctic Research Center
930 Koyukuk Drive, PO Box 757340
Fairbanks, Alaska 99775-7340
907-474-1939 office
907-385-7657 cell
kebenn...@alaska.edu
Personal Address:
UAF, PO Box 752525
Fairbanks, Alaska 99775-2525
bennett.katr...@gmail.com
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Timeseries Data Plotted as Monthly Boxplots
Hello, I'm trying to develop a box plot of time series data to look at the
range in the data values over the entire period of record.
My data initially starts out as a list of hourly data, and then I've been
using this code to make this data into the final ts array.
# Read in the station list
stn.list <- read.csv("/home/kbennett/fews/stnlist3", as.is=T, header=F)
# Read in all three variables.
vars <- c("MAT", "MAP", "MAP06")
for (stn in stn.list) {
for (v in 1:length(vars) {
# Read in year month start and end dates table & name it
ym.table <- read.csv("/home/kbennett/fews/", stn, var, ".ym.txt", as.is=T,
header=F)
names(ym.table) <- c("yearstart", "monthstart", "yearend", "monthend")
fn <- paste(stn, ".", vars[v], ".FIN", sep="")
if(file.exists(fn)) {
clim.dat <- read.csv(fn, header=F)
names(clim.dat) <- c("cdata")
year.start <- ym.table$yearstart
year.end <- ym.table$yearend
mo.start <- ym.table$monthstart
mo.end <- ym.table$monthend
regts.start = ISOdatetime(year.start, mo.start, 1, hour=0, min=0,
sec=0, tz="GMT")
regts.end = ISOdatetime(year.end, mo.end, 1, hour=18, min=0, sec=0,
tz="GMT")
zts <- zooreg(clim.dat$cdata, start = regts.start, end = regts.end,
frequency = 4, deltat = 21600)
#Create a daily average from the timeseries
zta <- aggregate(zts, as.POSIXct(cut(time(zts), "24 hours",
include=T)), mean)
#Select hourly data from the timeseries based on a specific time
zt.hr <- aggregate(zts, as.Date, head, 4)
zt.hr.ym <- aggregate(zt.hr, as.yearmon, head, 4)
zt.hr.1 <- zt.hr.ym[,1]
zt.hr.2 <- zt.hr.ym[,2]
zt.hr.3 <- zt.hr.ym[,3]
zt.hr.4 <- zt.hr.ym[,4]
zt.hr.1a <- aggregate(zt.hr.1, as.yearmon)
min.y <- min(zt.hr)
max.y <- max(zt.hr)
frequency(zt.hr.1) <- 12
zt.1.mo <- as.ts(zt.hr.1)
#Monthly boxplots of daily averages, for the months
boxplot(zt.1.mo ~ month, ##THIS IS WHAT DOESN'T
WORK HERE
boxwex=0.25, at=(1:12)-0.2,
outline = F,
col = "gray",
xlab = "Month",
ylab = expression(paste("( ",T^o,"C )") ),
ylim = c(min.y-5,max.y+5),
yaxs = "i",
xaxt = "n",
main = vars)
axis(1, at=c(1:12), labels=month.abb, cex.axis = 0.65)
legend("topright", c("Hour 00"), fill = c("gray"))
}
#write the results to a csv file
write.csv(cdat, paste(stn, "_", vars[v], ".csv", sep=""),
row.names=T, col.names=T)
}
}
The final array looks like this:
JanFebMarAprMayJunJulAugSep
OctNovDec
1948 28.719 4.977 39.037 9.746 8.348 36.672 47.660 54.076 38.062 34.486
11.938 39.666
1949 11.698 -6.675 16.844 0.950 10.349 38.752 39.785 40.544 57.603 35.476
2.308 -7.960
1950 0.340 45.206 6.385 17.132 19.074 38.465 48.711 54.686 48.743 33.978
23.090 10.007
1951 12.398 31.304 47.182 4.539 23.223 45.668 50.516 53.239 59.402 28.081
16.427 14.839
1952 -7.693 30.561 33.478 14.799 12.750 35.359 43.180 57.840 44.593 43.768
8.574 14.587
1953 -9.875 38.726 26.393 12.881 19.228 48.833 49.903 56.224 48.829 23.783
19.308 14.292
1954 35.943 16.706 16.021 7.806 23.593 40.418 45.310 53.113 49.203 29.480
17.228 33.068
1955 23.363 15.706 14.100 17.271 19.258 36.969 47.301 51.826 40.446 35.201
16.463 11.132
1956 45.868 -8.504 48.167 10.746 25.024 36.247 47.741 52.160 41.781 29.115
25.414 21.954
My main problem is that I can't access the rows (i.e. months) to subset the
data by.
Could someone point out how I am able to get at the months in this array and
subset them for plotting using the boxplot function?
Thank you,
Katrina
--
Katrina E. Bennett
PhD Student
University of Alaska Fairbanks
International Arctic Research Center
930 Koyukuk Drive, PO Box 757340
Fairbanks, Alaska 99775-7340
907-474-1939 office
907-385-7657 cell
kebenn...@alaska.edu
Personal Address:
UAF, PO Box 752525
Fairbanks, Alaska 99775-2525
bennett.katr...@gmail.com
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] zooreg and window
Hello, I have a following time series data
head(mend.dat)
ID PARAM Year Month Day Value SYM
1 15052500 1 1965 5 15 128 A
2 15052500 1 1965 5 16 135 A
3 15052500 1 1965 5 17 157 A
4 15052500 1 1965 5 18 176 A
5 15052500 1 1965 5 19 198 A
6 15052500 1 1965 5 20 241 A
I have loaded these data into as a zooreg object. The entire time
series runs from 1965-05-15 until 2010-12-31.
I have generated dates using ISOdatetime.
regts.start <- ISOdatetime(year.start, mo.start, day.start, hour=0,
min=0, sec=0, tz="GMT")
regts.end <- ISOdatetime(year.end, mo.end, day.end, hour=0, min=0,
sec=0, tz="GMT")
Then, I generate a zoo time series from the data.
zts <- zooreg(mend.dat$Value, start = regts.start, end = regts.end,
frequency = 1, deltat = 86400)
I now want to subset using window.
window(zts, start=as.Date("2000-01-01"), end=as.Date("2000-01-02"))
This results in the following error message.
Data:
factor(0)
1496 Levels: 100 1000 1 101 1010 10100 102 1020 103 1030 104 1040
105 1050 10500 106 1060 107 1070 108 1080 ... Eqp
Index:
character(0)
Warning messages:
1: In which(in.index & all.indexes >= start & all.indexes <= end) :
Incompatible methods ("Ops.POSIXt", "Ops.Date") for ">="
2: In which(in.index & all.indexes >= start & all.indexes <= end) :
Incompatible methods ("Ops.POSIXt", "Ops.Date") for "<="
Why is this happening?
Thanks for your help.
Katrina Bennett
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] zooreg and window
Hi Achim,
Yes, this worked. I added the two strings
regts.start <- as.Date(regts.start)
regts.end <- as.Date(regts.end)
and changed the frequency to 1 in my zooreg call.
The window call now works perfectly.
I'm still wrapping my head around POSIXct and ISOdatetime and how they
need to be applied appropriately!
Thank you,
Katrina
On Thu, Apr 14, 2011 at 12:26 AM, Achim Zeileis
wrote:
> On Wed, 13 Apr 2011, Katrina Bennett wrote:
>
>> Hello, I have a following time series data
>>
>> head(mend.dat)
>>
>> ID PARAM Year Month Day Value SYM
>> 1 15052500 1 1965 5 15 128 A
>> 2 15052500 1 1965 5 16 135 A
>> 3 15052500 1 1965 5 17 157 A
>> 4 15052500 1 1965 5 18 176 A
>> 5 15052500 1 1965 5 19 198 A
>> 6 15052500 1 1965 5 20 241 A
>>
>>
>> I have loaded these data into as a zooreg object. The entire time
>> series runs from 1965-05-15 until 2010-12-31.
>>
>> I have generated dates using ISOdatetime.
>>
>> regts.start <- ISOdatetime(year.start, mo.start, day.start, hour=0,
>> min=0, sec=0, tz="GMT")
>> regts.end <- ISOdatetime(year.end, mo.end, day.end, hour=0, min=0,
>> sec=0, tz="GMT")
>
> Here, you use "POSIXct" for your time index, below in the window() call you
> use "Date". The latter is more appropriate here because you appear to have
> daily data (and not intraday data).
>
> Thus, you can do
>
> regts.start <- as.Date(regts.start)
>
> or create it from scratch via
>
> regts.start <- as.Date("1965-05-15")
>
> and analogously for regts.end.
>
> The deltat in the zooreg() call is then simply 1.
>
>> Then, I generate a zoo time series from the data.
>>
>> zts <- zooreg(mend.dat$Value, start = regts.start, end = regts.end,
>> frequency = 1, deltat = 86400)
>>
>> I now want to subset using window.
>>
>> window(zts, start=as.Date("2000-01-01"), end=as.Date("2000-01-02"))
>>
>> This results in the following error message.
>>
>>
>> Data:
>> factor(0)
>> 1496 Levels: 100 1000 1 101 1010 10100 102 1020 103 1030 104 1040
>> 105 1050 10500 106 1060 107 1070 108 1080 ... Eqp
>>
>> Index:
>> character(0)
>> Warning messages:
>> 1: In which(in.index & all.indexes >= start & all.indexes <= end) :
>> Incompatible methods ("Ops.POSIXt", "Ops.Date") for ">="
>> 2: In which(in.index & all.indexes >= start & all.indexes <= end) :
>> Incompatible methods ("Ops.POSIXt", "Ops.Date") for "<="
>>
>>
>> Why is this happening?
>>
>> Thanks for your help.
>>
>> Katrina Bennett
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
--
Katrina E. Bennett
PhD Student
University of Alaska Fairbanks
International Arctic Research Center
930 Koyukuk Drive, PO Box 757340
Fairbanks, Alaska 99775-7340
907-474-1939 office
907-385-7657 cell
kebenn...@alaska.edu
Personal Address:
UAF, PO Box 752525
Fairbanks, Alaska 99775-2525
bennett.katr...@gmail.com
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Translate Sine Function in R?
Hello, I'm trying to generate a sine wave in R to fit my observations using the general formula: y=a*sin(b[x+h*pi)]+k where a = amplitude, b=period, h=phase shift, and k=vertical shift I want to use following translation to bring the sine function up onto the y-axis to range from 0-1, and this will place the wave on the x-axis from 0-pi/2. y=1/2sin(2[x+ 1/4*pi]) + 1/2 Additionally, I need to spread this along a x-axis that spans 1-153 (days). Can anyone help with this? I seem to be able to use the curve function fine, but entering the translations doesn't seem to provide an answer. Here is an example of the data set I am trying to 'match' using this function. dat <- c(75.44855206,NA,NA,NA,82.70745342,82.5335019,88.56617647,80.00128866,94.15418227,86.63987539,93.91052952,74.10612245,86.62289562,90.47961047,NA,NA,82.45320197,72.14371257,NA,71.44104803,72.59742896,68.36363636,NA,NA,61,NA,NA,71.26502909,NA,85.9333,84.34248284,79.00522193,79.64223058,97.2074017,88.43700548,96.40413877,95.13511869,92.57379057,93.97498475,NA,97.55995131,89.53321146,97.21728545,93.21980198,77.54054054,95.85392575,86.25684723,97.55325624,80.03950617,NA,91.34023128,92.42906574,88.59433962,65.77272727,89.63772455,NA,NA,NA,NA,74.86344239,83.57594937,70.22516556,65.30543319,NA,NA,67.84852294,60.90909091,54.79303797,NA,52.18735363,33.47003155,NA,41.34693878,24.5047043,NA,NA,NA,NA,9.9,13.6875,NA,11.90267176,84.14285714,3.781456954,NA,1.432926829,4.26557377,1.823529412,0.444620253,4.711155378,NA,6.320284698,0.581632653,0.144578313,3.7,0,0,0,0,0,NA,0.032947462,0,0,10.54545455,0,NA,0.561007958,0.75,NA,0.048780488,0.74137931,NA,2.023339318,0,0,0,NA,NA,0.156950673,NA,0.283769634,32.81818182,NA,NA,0,NA,0,0,0,NA,0.212454212,3.120181406,NA,0.011811024,NA,0,0.120430108,5.928571429,1.75,0.679292929,0.97,NA,0,NA,NA,1,0.38547486,NA,1.460732984,0.007795889,0.05465288,0.004341534) plot(dat/100) par(new=F) x.seq <- seq(100, 0, , 153) y <- ??? #*y = 2 sin 2Ï (x - 1/4)* or y ~ a + c*sin(x+b) However, I can't find a reference for the no place for k. Also, I've tried a lot of different iterations, but can't seem to figure out how to do this in R. Any thoughts or ideas on this? Thank you, Katrina [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Translate Sine Function in R?
Hi Petr, Thank you for the response. I am not familiar with the nls package. I am interested in obtaining a function to represent the relationship observed in "dat" over x.seq (days) so that I can calculate a minimum and maximum of the function. In the example below, it seems like you are using a data frame, but I tried to create the data frame from x.seq and dat: dat.df <- data.frame(x.seq, dat) However, I got the following error: Error in nls(dat ~ SSlogis(x, Asym, xmid, scal), data = dat.df, start = list(Asym = 90, : parameters without starting value in 'data': x Thank you for your assistance, Katrina On Thu, Aug 4, 2011 at 11:14 PM, Petr PIKAL wrote: > Hi > > Are you sure about sine fit? Seems to me that logistics would be better > > fit<-nls(dat ~ SSlogis(x, Asym, xmid,scal), data = dat.df, start = > list(Asym=90, xmid = 75, scal = -6)) > plot(dat.df) > lines(dat.df$x[complete.cases(dat.df)], predict(fit)) > > Regards > Petr > > > > > > > Hello, I'm trying to generate a sine wave in R to fit my observations > using > > the general formula: > > > > y=a*sin(b[x+h*pi)]+k > > > > where a = amplitude, b=period, h=phase shift, and k=vertical shift > > > > I want to use following translation to bring the sine function up onto > the > > y-axis to range from 0-1, and this will place the wave on the x-axis > from > > 0-pi/2. > > > > y=1/2sin(2[x+ 1/4*pi]) + 1/2 > > > > Additionally, I need to spread this along a x-axis that spans 1-153 > (days). > > > > Can anyone help with this? I seem to be able to use the curve function > fine, > > but entering the translations doesn't seem to provide an answer. > > > > Here is an example of the data set I am trying to 'match' using this > > function. > > > > dat <- > > > c(75.44855206,NA,NA,NA,82.70745342,82.5335019,88.56617647,80.00128866,94. > > 15418227,86.63987539,93.91052952,74.10612245,86.62289562,90. > > 47961047,NA,NA,82.45320197,72.14371257,NA,71.44104803,72.59742896,68. > > 36363636,NA,NA,61,NA,NA,71.26502909,NA,85.9333,84.34248284,79. > > 00522193,79.64223058,97.2074017,88.43700548,96.40413877,95.13511869,92. > > 57379057,93.97498475,NA,97.55995131,89.53321146,97.21728545,93.21980198, > > 77.54054054,95.85392575,86.25684723,97.55325624,80.03950617,NA,91. > > 34023128,92.42906574,88.59433962,65.77272727,89.63772455,NA,NA,NA,NA,74. > > 86344239,83.57594937,70.22516556,65.30543319,NA,NA,67.84852294,60. > > 90909091,54.79303797,NA,52.18735363,33.47003155,NA,41.34693878,24. > > 5047043,NA,NA,NA,NA,9.9,13.6875,NA,11.90267176,84.14285714,3. > > 781456954,NA,1.432926829,4.26557377,1.823529412,0.444620253,4. > > > 711155378,NA,6.320284698,0.581632653,0.144578313,3.7,0,0,0,0,0,NA, > > 0.032947462,0,0,10.54545455,0,NA,0.561007958,0.75,NA,0.048780488,0. > > 74137931,NA,2.023339318,0,0,0,NA,NA,0.156950673,NA,0.283769634,32. > > > 81818182,NA,NA,0,NA,0,0,0,NA,0.212454212,3.120181406,NA,0.011811024,NA,0, > > > 0.120430108,5.928571429,1.75,0.679292929,0.97,NA,0,NA,NA,1,0.38547486,NA, > > 1.460732984,0.007795889,0.05465288,0.004341534) > > plot(dat/100) > > par(new=F) > > x.seq <- seq(100, 0, , 153) > > y <- ??? #*y = 2 sin 2Ä⬠(x - 1/4)* or y ~ a + c*sin(x+b) > > > > However, I can't find a reference for the no place for k. Also, I've > tried a > > lot of different iterations, but can't seem to figure out how to do this > in > > R. > > > > Any thoughts or ideas on this? > > > > Thank you, > > > > Katrina > > > >[[alternative HTML version deleted]] > > > > __ > > R-help@r-project.org mailing list > > https://stat.ethz.ch/mailman/listinfo/r-help > > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > > and provide commented, minimal, self-contained, reproducible code. > > -- Katrina E. Bennett PhD Student University of Alaska Fairbanks International Arctic Research Center 930 Koyukuk Drive, PO Box 757340 Fairbanks, Alaska 99775-7340 907-474-1939 office 907-385-7657 cell kebenn...@alaska.edu Personal Address: UAF, PO Box 752525 Fairbanks, Alaska 99775-2525 bennett.katr...@gmail.com [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] Translate Sine Function in R?
Sorry, I sent that e-mail too soon. I realize you called "x.seq" just "x" in the data frame. Yes, this is a great way to fit these data. I think from this that I can probably now work with the related function of the line to get my minimum and maximum values. Thanks for your assistance! Katrina On Fri, Aug 5, 2011 at 4:41 PM, Katrina Bennett wrote: > Hi Petr, > Thank you for the response. > > I am not familiar with the nls package. I am interested in obtaining a > function to represent the relationship observed in "dat" over x.seq (days) > so that I can calculate a minimum and maximum of the function. > > In the example below, it seems like you are using a data frame, but I tried > to create the data frame from x.seq and dat: > > dat.df <- data.frame(x.seq, dat) > > However, I got the following error: > > Error in nls(dat ~ SSlogis(x, Asym, xmid, scal), data = dat.df, start = > list(Asym = 90, : > parameters without starting value in 'data': x > > Thank you for your assistance, > > Katrina > > > On Thu, Aug 4, 2011 at 11:14 PM, Petr PIKAL wrote: > >> Hi >> >> Are you sure about sine fit? Seems to me that logistics would be better >> >> fit<-nls(dat ~ SSlogis(x, Asym, xmid,scal), data = dat.df, start = >> list(Asym=90, xmid = 75, scal = -6)) >> plot(dat.df) >> lines(dat.df$x[complete.cases(dat.df)], predict(fit)) >> >> Regards >> Petr >> >> >> >> > >> > Hello, I'm trying to generate a sine wave in R to fit my observations >> using >> > the general formula: >> > >> > y=a*sin(b[x+h*pi)]+k >> > >> > where a = amplitude, b=period, h=phase shift, and k=vertical shift >> > >> > I want to use following translation to bring the sine function up onto >> the >> > y-axis to range from 0-1, and this will place the wave on the x-axis >> from >> > 0-pi/2. >> > >> > y=1/2sin(2[x+ 1/4*pi]) + 1/2 >> > >> > Additionally, I need to spread this along a x-axis that spans 1-153 >> (days). >> > >> > Can anyone help with this? I seem to be able to use the curve function >> fine, >> > but entering the translations doesn't seem to provide an answer. >> > >> > Here is an example of the data set I am trying to 'match' using this >> > function. >> > >> > dat <- >> > >> c(75.44855206,NA,NA,NA,82.70745342,82.5335019,88.56617647,80.00128866,94. >> > 15418227,86.63987539,93.91052952,74.10612245,86.62289562,90. >> > 47961047,NA,NA,82.45320197,72.14371257,NA,71.44104803,72.59742896,68. >> > 36363636,NA,NA,61,NA,NA,71.26502909,NA,85.9333,84.34248284,79. >> > 00522193,79.64223058,97.2074017,88.43700548,96.40413877,95.13511869,92. >> > 57379057,93.97498475,NA,97.55995131,89 >> .53321146,97.21728545,93.21980198, >> > 77.54054054,95.85392575,86.25684723,97.55325624,80.03950617,NA,91. >> > 34023128,92.42906574,88.59433962,65.77272727,89.63772455,NA,NA,NA,NA,74. >> > 86344239,83.57594937,70.22516556,65.30543319,NA,NA,67.84852294,60. >> > 90909091,54.79303797,NA,52.18735363,33.47003155,NA,41.34693878,24. >> > 5047043,NA,NA,NA,NA,9.9,13.6875,NA,11.90267176,84.14285714,3. >> > 781456954,NA,1.432926829,4.26557377,1.823529412,0.444620253,4. >> > >> 711155378,NA,6.320284698,0.581632653,0.144578313,3.7,0,0,0,0,0,NA, >> > 0.032947462,0,0,10.54545455,0,NA,0.561007958,0.75,NA,0.048780488,0. >> > 74137931,NA,2.023339318,0,0,0,NA,NA,0.156950673,NA,0.283769634,32. >> > >> 81818182,NA,NA,0,NA,0,0,0,NA,0.212454212,3.120181406,NA,0.011811024,NA,0, >> > >> 0.120430108,5.928571429,1.75,0.679292929,0.97,NA,0,NA,NA,1,0.38547486,NA, >> > 1.460732984,0.007795889,0.05465288,0.004341534) >> > plot(dat/100) >> > par(new=F) >> > x.seq <- seq(100, 0, , 153) >> > y <- ??? #*y = 2 sin 2Ä⬠(x - 1/4)* or y ~ a + c*sin(x+b) >> > >> > However, I can't find a reference for the no place for k. Also, I've >> tried a >> > lot of different iterations, but can't seem to figure out how to do this >> in >> > R. >> > >> > Any thoughts or ideas on this? >> > >> > Thank you, >> > >> > Katrina >> > >> >[[alternative HTML version deleted]] >> > >> > __ >> > R-help@r-project.org mailing list >> > https://stat.ethz.ch/mailman/listinfo/r-help >> > PLEASE do read the posting guide >> http://www.R-project.org/posting-guide.html >> > and provide commented, minimal, self-contained, reproducible code. >> >> [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] nls, how to determine function?
Hi R help, I am trying to determine how nls() generates a function based on the self-starting SSlogis and what the formula for the function would be. I've scoured the help site, and other literature to try and figure this out but I still am unsure if I am correct in what I am coming up with. ** dat <- c(75.44855206,NA,NA,NA,82.70745342,82.5335019,88.56617647,80.00128866,94.15418227,86.63987539,93.91052952,74.10612245,86.62289562,90.47961047,NA,NA,82.45320197,72.14371257,NA,71.44104803,72.59742896,68.36363636,NA,NA,61,NA,NA,71.26502909,NA,85.9333,84.34248284,79.00522193,79.64223058,97.2074017,88.43700548,96.40413877,95.13511869,92.57379057,93.97498475,NA,97.55995131,89.53321146,97.21728545,93.21980198,77.54054054,95.85392575,86.25684723,97.55325624,80.03950617,NA,91.34023128,92.42906574,88.59433962,65.77272727,89.63772455,NA,NA,NA,NA,74.86344239,83.57594937,70.22516556,65.30543319,NA,NA,67.84852294,60.90909091,54.79303797,NA,52.18735363,33.47003155,NA,41.34693878,24.5047043,NA,NA,NA,NA,9.9,13.6875,NA,11.90267176,84.14285714,3.781456954,NA,1.432926829,4.26557377,1.823529412,0.444620253,4.711155378,NA,6.320284698,0.581632653,0.144578313,3.7,0,0,0,0,0,NA,0.032947462,0,0,10.54545455,0,NA,0.561007958,0.75,NA,0.048780488,0.74137931,NA,2.023339318,0,0,0,NA,NA,0.156950673,NA,0.283769634,32.81818182,NA,NA,0,NA,0,0,0,NA,0.212454212,3.120181406,NA,0.011811024,NA,0,0.120430108,5.928571429,1.75,0.679292929,0.97,NA,0,NA,NA,1,0.38547486,NA,1.460732984,0.007795889,0.05465288,0.004341534) dat.df.1 <- data.frame(dat) dat.df.2 <- data.frame(x=x.seq, dat.df=dat.df.1) fit.dat <-nls(dat ~ SSlogis(x, Asym, xmid,scal), data = dat.df.2, start =list(Asym=90, xmid = 75, scal = -6)) plot(dat.df.2, axes=FALSE, ann=FALSE, ylim=c(0,100)) lines(dat.df.2$x[complete.cases(dat.df.2)], predict(fit.dat), ylim=c(0,100)) summary(fit.dat) ** Formula: dat ~ SSlogis(x, Asym, xmid, scal) Parameters: Estimate Std. Error t value Pr(>|t|) Asym 85.651 1.716 49.900 < 2e-16 *** xmid 72.214 1.036 69.697 < 2e-16 *** scal -6.150 0.850 -7.236 7.9e-11 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 10.33 on 105 degrees of freedom Number of iterations to convergence: 10 Achieved convergence tolerance: 4.405e-06 (45 observations deleted due to missingness) ** >From r-help, SSlogis parameters asym, xmid and scal are defined as: Asym: a numeric parameter representing the asymptote. xmid: a numeric parameter representing the x value at the inflection point of the curve. The value of SSlogis will be Asym/2 at xmid. scal: a numeric scale parameter on the input axis. and it states that the value of SSlogis "is a numeric vector of the same length as input. It is the value of the expression sym/(1+exp((xmid-input)/scal)). If all of the arguments Asym, xmid, and scal are names of objects the gradient matrix with respect to these names is attached as an attribute named gradient." However, how do I get the actual function for the curve that is generated? I don't think it can just be: y= asym/((1+e^((xmid-x)/scal)))? Also, how do you determine the starting parameters to input in for asym, xmin, and scal? Perhaps I need to start at the beginning and define my own function, and not rely on SSlogis to provide it? What I want to be able to do is determine a local maximum for my curve (the x value at which this curve inflects (the upper inflection)), and the x value for the local minimum (the lower inflection curve), and the x value counts in between these values. I think in order to do this I need to differentiate the function. Any insight on this would be greatly appreciated. Sincerely, Katrina __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] nls, how to determine function?
Hi Petr, thanks for your help on this. I will most definitely get this book as it appears to be a good one. I am happy with how nls() and the self starting function appears to be fitting the data set. What I am wondering about is how to write out this function outside of R. For example, if I sat down with a paper and pencil and wrote this function out in terms of f(x) = ?. In this example (http://www.graphpad.com/curvefit/id203.htm) there is a formula contains for what they are calling the Boltzmann Sigmoid function. In order to have an answer, is this something I need to generate on my own? Or, is there a common form that is being applied by the SSlogis and the SSfpl functions? In Pinheiro and Bates (pg 518) there is a model f(x) = A + B/1+ exp[(scal-x)/expL]. Is this the form of the equation? Or, is there a way to extract the local min and max values from nls () results? Do I have to break apart the individual components of the regression into piecewise functions? Thanks again for your help. Katrina On Tue, Aug 9, 2011 at 12:38 AM, Petr PIKAL wrote: > Hi > >> Hi R help, >> >> I am trying to determine how nls() generates a function based on the >> self-starting SSlogis and what the formula for the function would be. >> I've scoured the help site, and other literature to try and figure >> this out but I still am unsure if I am correct in what I am coming up >> with. > > Thanks for providing data and your code > >> >> >> > ** >> dat <- c(75.44855206,NA,NA,NA,82.70745342,82.5335019,88.56617647,80. >> 00128866,94.15418227,86.63987539,93.91052952,74.10612245,86.62289562,90. >> 47961047,NA,NA,82.45320197,72.14371257,NA,71.44104803,72.59742896,68. >> 36363636,NA,NA,61,NA,NA,71.26502909,NA,85.9333,84.34248284,79. >> 00522193,79.64223058,97.2074017,88.43700548,96.40413877,95.13511869,92. >> 57379057,93.97498475,NA,97.55995131,89.53321146,97.21728545,93.21980198, >> 77.54054054,95.85392575,86.25684723,97.55325624,80.03950617,NA,91. >> 34023128,92.42906574,88.59433962,65.77272727,89.63772455,NA,NA,NA,NA,74. >> 86344239,83.57594937,70.22516556,65.30543319,NA,NA,67.84852294,60. >> 90909091,54.79303797,NA,52.18735363,33.47003155,NA,41.34693878,24. >> 5047043,NA,NA,NA,NA,9.9,13.6875,NA,11.90267176,84.14285714,3. >> 781456954,NA,1.432926829,4.26557377,1.823529412,0.444620253,4. >> > 711155378,NA,6.320284698,0.581632653,0.144578313,3.7,0,0,0,0,0,NA, >> 0.032947462,0,0,10.54545455,0,NA,0.561007958,0.75,NA,0.048780488,0. >> 74137931,NA,2.023339318,0,0,0,NA,NA,0.156950673,NA,0.283769634,32. >> > 81818182,NA,NA,0,NA,0,0,0,NA,0.212454212,3.120181406,NA,0.011811024,NA,0, >> > 0.120430108,5.928571429,1.75,0.679292929,0.97,NA,0,NA,NA,1,0.38547486,NA, >> 1.460732984,0.007795889,0.05465288,0.004341534) > >> dat.df.1 <- data.frame(dat) > unnecessary > >> dat.df.2 <- data.frame(x=x.seq, dat.df=dat.df.1) > > some correction > dat.df.2 <- data.frame(x=seq_along(dat), dat=dat) > >> fit.dat <-nls(dat ~ SSlogis(x, Asym, xmid,scal), data = dat.df.2, >> start =list(Asym=90, xmid = 75, scal = -6)) >> plot(dat.df.2, axes=FALSE, ann=FALSE, ylim=c(0,100)) >> lines(dat.df.2$x[complete.cases(dat.df.2)], predict(fit.dat), > ylim=c(0,100)) >> >> summary(fit.dat) >> >> > ** >> Formula: dat ~ SSlogis(x, Asym, xmid, scal) >> >> Parameters: >> Estimate Std. Error t value Pr(>|t|) >> Asym 85.651 1.716 49.900 < 2e-16 *** >> xmid 72.214 1.036 69.697 < 2e-16 *** >> scal -6.150 0.850 -7.236 7.9e-11 *** >> --- >> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 >> >> Residual standard error: 10.33 on 105 degrees of freedom >> >> Number of iterations to convergence: 10 >> Achieved convergence tolerance: 4.405e-06 >> (45 observations deleted due to missingness) >> > ** >> >> >From r-help, SSlogis parameters asym, xmid and scal are defined as: >> >> Asym: a numeric parameter representing the asymptote. >> >> xmid: a numeric parameter representing the x value at the inflection >> point of the curve. The value of SSlogis will be Asym/2 at xmid. >> >> scal: a numeric scale parameter on the input axis. >> >> and it states that the value of SSlogis "is a numeric vector of the >> same length as input. It is the value of the expression >> sym/(1+exp((xmid-input)/scal)). If all of the arguments Asym, xmid, >> and scal are names of objects the gradient matrix with respect to >> these names is attached as an attribute named gradient." >> >> However, how do I get the actual function for the curve that is >> generated? I don't think it can just be: y= >> asym/((1+e^((xmid-x)/scal)))? > > Yes. I think that the best source of information about nonlinear > regression is book by Bates, Pinheiro - Mixed effect models with S and S+. > There you can find how to determine starting parameters, h
[R] sapply to bind columns, with repeat?
Hi R-help,
I am working with US COOP network station data and the files are
concatenated in single rows for all years, but I need to pull these
apart into rows for each day. To do this, I need to extract part of
each row such as station id, year, mo, and repeat this against other
variables in the row (days). My problem is that there are repeated
values for each day, and the files are fixed width field without
order.
Here is an example of just one line of data.
coop.raw <- c("DLY09752806TMAX F2010010620107 00049 20107 00062
B0207 00041 20207 00049 B0307 00040 20307 00041 B0407 00042 20407
00040 B0507 00041 20507 00042 B0607 00043 20607 00041 B0707 00055
20707 00043 B0807 00039 20807 00055 B0907 00037 20907 00039 B1007
00038 21007 00037 B1107 00048 21107 00038 B1207 00050 21207 00048
B1307 00051 21307 00050 B1407 00058 21407 00051 B1507 00068 21507
00058 B1607 00065 21607 00068 B1707 00068 21707 00065 B1807 00067
21807 00068 B1907 00068 21907 00067 B2007 00069 22007 00068 B2107
00057 22107 00069 B2207 00048 22207 00057 B2307 00051 22307 00048
B2407 00073 22407 00051 B2507 00062 22507 00073 B2607 00056 22607
00062 B2707 00053 22707 00056 B2807 00064 22807 00053 B2907 00057
22907 00064 B3007 00047 23007 00057 B3107 00046 23107 00047 B")
write.csv(coop.raw, "coop.tmp", row.names=F, quote=F)
coop.dat <- read.fwf("coop.tmp", widths =
c(c(3,8,4,2,4,2,4,3),rep(c(2,2,1,5,1,1),62)), na.strings=c(""),
skip=1, as.is=T)
rep.name <- rep(c("day","hr","met","dat","fl1","fl2"), 62)
rep.count <- rep(c(1:62), each=6, 1)
names(coop.dat) <- c("rect", "id", "elem", "unt", "year", "mo",
"fill", "numval", paste(rep.name, rep.count, sep="_"))
I would like to generate output that contains in one row, the columns
"id", "elem", "unt", "year", "mo", and "numval". Binded to these
initial columns, I would like only "day_1", "hr_1", "met_1", "dat_1",
"fl1_1", and "fl2_1". Then, in the next row I would like repeated the
initial columns "id", "elem", "unt", "year", "mo", and "numval" and
then binded "day_2", "hr_2", "met_2", "dat_2", "fl1_2", and "f2_2" and
so on until all the data for all rows has been allocated. Then, move
onto the next row and repeat.
I think I should be able to do this with some sort of sapply or lapply
function, but I'm struggling with the format for repeating the initial
columns, and then skipping through the next columns.
Thank you,
Katrina
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] sapply to bind columns, with repeat?
Hi Weidong Gu,
This works! For my clarity, and so I can repeat this process if need be:
The 'mat' generates a matrix using whatever is supplied to x (i.e.
coop.dat) using the columns from position 9:length(x) of 6 columns (by
row).
The 'rem.col' generates a matrix of the first 1:8 columns of 8 columns.
The 'return' statement calls the function to cbind together rem.col and mat.
Then 'apply' this all to coop.dat, by rows, using function reorg.
Is this correct?
Thank you very much,
Katrina
On Fri, Aug 12, 2011 at 10:28 AM, Weidong Gu wrote:
> Katrina,
>
> try this.
>
> reorg<-function(x){
> mat<-matrix(x[9:length(x)],ncol=6,byrow=T)
> rem.col<-matrix(rep(x[1:8],nrow(mat)),byrow=T,ncol=8)
> return(data.frame(cbind(rem.col,mat)))
> }
>
> co<-do.call('rbind',apply(coop.dat,1,function(x) reorg(x)))
>
> You may need to tweak a bit to fit exactly what you want.
>
> Weidong Gu
>
> On Fri, Aug 12, 2011 at 2:35 AM, Katrina Bennett wrote:
>> Hi R-help,
>>
>> I am working with US COOP network station data and the files are
>> concatenated in single rows for all years, but I need to pull these
>> apart into rows for each day. To do this, I need to extract part of
>> each row such as station id, year, mo, and repeat this against other
>> variables in the row (days). My problem is that there are repeated
>> values for each day, and the files are fixed width field without
>> order.
>>
>> Here is an example of just one line of data.
>>
>> coop.raw <- c("DLY09752806TMAX F2010010620107 00049 20107 00062
>> B0207 00041 20207 00049 B0307 00040 20307 00041 B0407 00042 20407
>> 00040 B0507 00041 20507 00042 B0607 00043 20607 00041 B0707 00055
>> 20707 00043 B0807 00039 20807 00055 B0907 00037 20907 00039 B1007
>> 00038 21007 00037 B1107 00048 21107 00038 B1207 00050 21207 00048
>> B1307 00051 21307 00050 B1407 00058 21407 00051 B1507 00068 21507
>> 00058 B1607 00065 21607 00068 B1707 00068 21707 00065 B1807 00067
>> 21807 00068 B1907 00068 21907 00067 B2007 00069 22007 00068 B2107
>> 00057 22107 00069 B2207 00048 22207 00057 B2307 00051 22307 00048
>> B2407 00073 22407 00051 B2507 00062 22507 00073 B2607 00056 22607
>> 00062 B2707 00053 22707 00056 B2807 00064 22807 00053 B2907 00057
>> 22907 00064 B3007 00047 23007 00057 B3107 00046 23107 00047 B")
>> write.csv(coop.raw, "coop.tmp", row.names=F, quote=F)
>> coop.dat <- read.fwf("coop.tmp", widths =
>> c(c(3,8,4,2,4,2,4,3),rep(c(2,2,1,5,1,1),62)), na.strings=c(""),
>> skip=1, as.is=T)
>> rep.name <- rep(c("day","hr","met","dat","fl1","fl2"), 62)
>> rep.count <- rep(c(1:62), each=6, 1)
>> names(coop.dat) <- c("rect", "id", "elem", "unt", "year", "mo",
>> "fill", "numval", paste(rep.name, rep.count, sep="_"))
>>
>> I would like to generate output that contains in one row, the columns
>> "id", "elem", "unt", "year", "mo", and "numval". Binded to these
>> initial columns, I would like only "day_1", "hr_1", "met_1", "dat_1",
>> "fl1_1", and "fl2_1". Then, in the next row I would like repeated the
>> initial columns "id", "elem", "unt", "year", "mo", and "numval" and
>> then binded "day_2", "hr_2", "met_2", "dat_2", "fl1_2", and "f2_2" and
>> so on until all the data for all rows has been allocated. Then, move
>> onto the next row and repeat.
>>
>> I think I should be able to do this with some sort of sapply or lapply
>> function, but I'm struggling with the format for repeating the initial
>> columns, and then skipping through the next columns.
>>
>> Thank you,
>>
>> Katrina
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Seasonal and 11-day subset for zoo object
I have a zooreg object and I want to be able to generate a value for seasons
and 11-day composites paste it onto my zoo data frame, along with year,
month and days.
Right now I have the following to work from:
eg. dat.zoo.mdy <- with(month.day.year(time(dat.zoo)), cbind(dat.zoo, year,
month, day, quarter = (month - 1) %/% 3 + 1, dow =
as.numeric(format(time(dat.zoo), "%w"
For the seasons, I have been trying to replace 'quarter' with a seasonal
value of "1" for Dec-Jan-Feb, "2" for Mar-Apr-May, "3" for Jun-Jul-Aug, "4"
for Sep-Oct-Nov.
dat.zoo.mdy <- with(month.day.year(time(dat.zoo)), cbind(dat.zoo, year,
month, day,
season=for(i in nrow(dat.zoo.mdy)) {
if (month[i] == 12) {
quarter[i]=1
} else if (month[i] == 3) {
quarter[i]=2
} else if (month[i] == 6) {
quarter[i]=3
} else quarter[i]=4 }, dow = as.numeric(format(time(dat.zoo),
"%w"
However, this gives me the error: "Error in zoo(structure(x, dim = dim(x)),
index(x), ...) :
x : attempt to define illegal zoo object"
I'd like to get an 11-day value as well to replace the dow in the first
example, but I'm still trying to figure out if there is an easy way to do
this in zoo.
Any ideas would be greatly appreciated.
Katrina
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Obtaining a derivative of nls() SSlogis function
Hello, I am wondering if someone can help me. I have the following function
that I derived using nls() SSlogis. I would like to find its derivative. I
thought I had done this using deriv(), but for some reason this isn't
working out for me.
Here is the function:
asym <- 84.951
xmid <- 66.90742
scal <- -6.3
x.seq <- seq(1, 153,, 153)
nls.fn <- asym/((1+exp((xmid-x.seq)/scal)))
try #1
deriv(nls.fn)
#get an Error in .Internal(deriv.default(expr, namevec, function.arg, tag,
hessian)) : 'namevec' is missing
try #2
deriv(nls.fn, namevec=c("asym", "xmid", "scal"))
#this doesn't seem to give me the expression, and the gradients are zero.
I've tried to do this with Ryacas as well, but I'm lost.
Can anyone help?
Thank you,
Katrina
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Obtaining a derivative of nls() SSlogis function
Thank you to all who have made great suggestions. Both Gabor and David's
methods work well for me. I see now where I was going wrong with this, I
wasn't loading it to a function. I also did not know expression could be
used in that way.
Thank you very much for your help!
Katrina
On Fri, Nov 18, 2011 at 6:02 AM, Gabor Grothendieck wrote:
> On Thu, Nov 17, 2011 at 4:40 PM, Katrina Bennett
> wrote:
> > Hello, I am wondering if someone can help me. I have the following
> function
> > that I derived using nls() SSlogis. I would like to find its derivative.
> I
> > thought I had done this using deriv(), but for some reason this isn't
> > working out for me.
> >
> > Here is the function:
> > asym <- 84.951
> > xmid <- 66.90742
> > scal <- -6.3
> >
> > x.seq <- seq(1, 153,, 153)
> > nls.fn <- asym/((1+exp((xmid-x.seq)/scal)))
> >
> > try #1
> > deriv(nls.fn)
> > #get an Error in .Internal(deriv.default(expr, namevec, function.arg,
> tag,
> > hessian)) : 'namevec' is missing
> >
> > try #2
> > deriv(nls.fn, namevec=c("asym", "xmid", "scal"))
> > #this doesn't seem to give me the expression, and the gradients are zero.
> >
> > I've tried to do this with Ryacas as well, but I'm lost.
> >
> > Can anyone help?
> >
>
> You can do that with plain R:
>
> > e <- quote(asym/((1+exp((xmid-x.seq)/scal
> > for(v in all.vars(e)) cat("deriv wrt", v, "is", format(D(e, v)), "\n")
> deriv wrt asym is 1/((1 + exp((xmid - x.seq)/scal)))
> deriv wrt xmid is -(asym * (exp((xmid - x.seq)/scal) * (1/scal))/((1 +
> exp((xmid - x.seq)/scal)))^2)
> deriv wrt x.seq is asym * (exp((xmid - x.seq)/scal) * (1/scal))/((1 +
> exp((xmid - x.seq)/scal)))^2
> deriv wrt scal is asym * (exp((xmid - x.seq)/scal) * ((xmid -
> x.seq)/scal^2))/((1 + exp((xmid - x.seq)/scal)))^2
>
>
> --
> Statistics & Software Consulting
> GKX Group, GKX Associates Inc.
> tel: 1-877-GKX-GROUP
> email: ggrothendieck at gmail.com
>
--
Katrina E. Bennett
PhD Student
University of Alaska Fairbanks
International Arctic Research Center
930 Koyukuk Drive, PO Box 757340
Fairbanks, Alaska 99775-7340
907-474-1939 office
907-385-7657 cell
kebenn...@alaska.edu
Personal Address:
UAF, PO Box 752525
Fairbanks, Alaska 99775-2525
bennett.katr...@gmail.com
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] Retain parts of a matrix
Hi all,
I'm working to apply a function that will generate a matrix of results only
when a specific criteria is met.
I want my final results to be a matrix with both the values that meet the
criteria (the results of the function), and those that to do in the same
positions in the matrix (the original numbers).
Here's a sample of what I would like to do:
t.mean.1.c <- c(-15, -20, -30, -20, -25, -35, -40, -8, 9, 10)
cdem <- c(300, 400, 700, 900, 1000, 250, 200, 300, 500, 650, 650, 1200,
1400, 3200, 2000)
cdem.mat <- matrix(NA, 10, 15)
cdem.mat <- cdem.mat(sapply(cdem[which(cdem >= 900)], function(x)
get("t.mean.1.c") - (x * -0.0065)))
I want cdem.mat to be the same size as the the first call of it, with NA
values where the cdem is < 899. I thought this would fill it in but it
seems to just over write it and I only end up with the rows where cdem is
greater than 900.
Thank you,
Katrina
--
Katrina E. Bennett
PhD Student
University of Alaska Fairbanks
International Arctic Research Center
930 Koyukuk Drive, PO Box 757340
Fairbanks, Alaska 99775-7340
907-474-1939 office
907-385-7657 cell
kebenn...@alaska.edu
Personal Address:
UAF, PO Box 752525
Fairbanks, Alaska 99775-2525
bennett.katr...@gmail.com
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Retain parts of a matrix
Sorry for not being more clear.
I'll try to explain again.
I have a rather large DEM and I need to scale daily temperature values for
10 years.
I am using the sapply function to adjust temperatures (t.mean.1.c) based on
lapse rates across DEM points (cdem) which meet the condition of elevation
being greater than 300m. Below ~300m the lapse rate calculate changes
because temperatures decrease with elevation. Above ~300m the lapse rates
calculation must account for an inversion, where temperature increases with
elevation.
What I would like as a result is the values of the DEM which are lower than
300m to retain a generic lapse rate adjustment, and the DEM values great
than 300 to be treating using a different function.
t.mean.1.c <- c(-15, -20, -30, -20, -25, -35, -40, -8, 9, 10)
cdem <- c(300, 400, 700, 900, 1000, 250, 200, 300, 500, 650, 650, 1200,
1400, 3200, 2000)
cdem.mat <- sapply(cdem[which(cdem >= 900)], function(x) get("t.mean.1.c")
- (x * -0.0065)) #this just gives me a
#result of only the six values
I'm going to make this even more simple to show what I want as output.
t.mean.1.c <- c(-14, -20)
cdem <- c(300, 400, 700, 900)
cdem.mat.1 <- sapply(cdem[which(cdem >= 900)], function(x)
get("t.mean.1.c") - (x * -0.0065))
cdem.mat.2 <- sapply(cdem[which(cdem < 900)], function(x) get("t.mean.1.c")
+ (x * -0.0065))
cdem.mat.1 #right now I get this, which is only the two temperature scaled
to the 900 m elevation case
# [,1]
#[1,] -8.15
#[2,] -14.15
#What I want:
# [,1][,2] [,3] [,4]
#[1,] -15.95 -16.6-18.55 -8.15
#[2,] -21.95 -22.6 -24.55 -14.15
Maybe I can combine in the sapply function an if else statement? I haven't
tested the below try.
sapply( if (cdem < 900)] function(x) get("t.mean.1.c") - (x * -0.0065)
else function(x) get("t.mean.1.c") + (x * -0.0065) )
Finally, I want to output everything into one single matrix array.
I hope that is more clear.
Thanks for your assistance.
Katrina
On Sun, Nov 27, 2011 at 7:35 PM, David Winsemius wrote:
>
> On Nov 27, 2011, at 10:15 PM, Katrina Bennett wrote:
>
> Hi all,
>>
>> I'm working to apply a function that will generate a matrix of results
>> only
>> when a specific criteria is met.
>>
>> I want my final results to be a matrix with both the values that meet the
>> criteria (the results of the function), and those that to do in the same
>> positions in the matrix (the original numbers).
>>
>> Here's a sample of what I would like to do:
>>
>> t.mean.1.c <- c(-15, -20, -30, -20, -25, -35, -40, -8, 9, 10)
>> cdem <- c(300, 400, 700, 900, 1000, 250, 200, 300, 500, 650, 650, 1200,
>> 1400, 3200, 2000)
>> cdem.mat <- matrix(NA, 10, 15)
>> cdem.mat <- cdem.mat(sapply(cdem[which(**cdem >= 900)], function(x)
>>
> ..^
> cdem.mat is not a function, cannot use "(".
>
>
>
> get("t.mean.1.c") - (x * -0.0065)))
>>
>> I want cdem.mat to be the same size as the the first call of it, with NA
>> values where the cdem is < 899.
>>
>
> You are only offering 10 values as replacements and an unequal number of
> condition regardless of whether you count the qualifying cdem conditions
> (6) or the total conditions(15). So you will need to be more clear about
> how you want the conditions to apply to values. They stepping the audience
> through the first few decisions and say where you want the values placed.
> Maybe you should also say whether you expect values to be recycled (the
> default with matrix assignment).
>
> I suppose it is possible that you want the results of this:
>
> cdem.mat[cdem>899] <- t.mean.1.c - (cdem * -0.0065)
>
> but I'm guessing not.
>
> I thought this would fill it in but it
>> seems to just over write it and I only end up with the rows where cdem is
>> greater than 900.
>>
>> Thank you,
>>
>> Katrina
>>
>>
>>
>> --
>> Katrina E. Bennett
>> PhD Student
>> University of Alaska Fairbanks
>> International Arctic Research Center
>> 930 Koyukuk Drive, PO Box 757340
>> Fairbanks, Alaska 99775-7340
>> 907-474-1939 office
>> 907-385-7657 cell
>> kebenn...@alaska.edu
>>
>>
>> Personal Address:
>> UAF, PO Box 752525
>> Fairbanks, Alaska 99775-2525
>> bennett.katr...@gmail.com
>>
>>[[alternative HTML version deleted]]
>>
>> __**
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/**listinfo/r-help<https://stat.ethz.ch/mailman/listinfo/r-help>
>> PLEASE do read the posting guide http://www.R-project.org/**
>> posting-guide.html <http://www.R-project.org/posting-guide.html>
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
> David Winsemius, MD
> West Hartford, CT
>
>
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Retain parts of a matrix
Hi David and Jim! Thanks very much for your help, this worked for me. #generate two matrices that are the same size for elevation and temperature: t.mean.1.c <- matrix(c(-15, -20, -30, -20, -25, -35, -40, -8, 9, 10), nrow=10, ncol=15, byrow=F) cdem <- matrix(c(300, 400, 700, 900, 1000, 250, 200, 300, 500, 650, 650, 1200, 1400, 3200, 2000),nrow=nrow(t.mean.1.c), ncol=15, byrow=T) #create the functions to use: gt300 <- function(elevation,temp) return((elevation * 0.0065) + temp) lt300 <- function(elevation,temp) return((elevation * -0.0065) + temp) #apply the function and subset everything by the values I want to apply it to cdem[cdem>=300] <- gt300(cdem[cdem>=300],t.mean.1.c[cdem>=300]) cdem[cdem<300] <- lt300(cdem[cdem<300],t.mean.1.c[cdem<300]) This retains the matrix and scales everything properly. Thanks again! Katrina p.s. Jim - yes, the dummy subset value I was using 900m but I switched it in my example to my actual value (since I am going to use 300m for this initial attempt at this). I will need to make it a lot more complicated eventually! On Mon, Nov 28, 2011 at 6:57 PM, Katrina Bennett wrote: > Hi David and Jim! > > Thanks very much for your help, this worked for me. > > #generate two matrices that are the same size for elevation and > temperature: > t.mean.1.c <- matrix(c(-15, -20, -30, -20, -25, -35, -40, -8, 9, 10), > nrow=10, ncol=15, byrow=F) > cdem <- matrix(c(300, 400, 700, 900, 1000, 250, 200, 300, 500, 650, 650, > 1200, 1400, 3200, 2000),nrow=nrow(t.mean.1.c),ncol=15, byrow=T) > > #create the functions to use: > gt300 <- function(elevation,temp) return(temp + (elevation * 0.0065)) > lt300 <- function(elevation,temp) return(temp + (elevation * -0.0065)) > > #apply the function and subset everything by the values I want to apply it > to > cdem[cdem>=300] <- gt300(cdem[cdem>=300],t.mean.1.c[cdem>=300]) > cdem[cdem<300] <- lt300(cdem[cdem<300],t.mean.1.c[cdem<300]) > > This retains the matrix and scales everything properly. > > Thanks again! > > Katrina > > p.s. Jim - yes, the dummy subset value I was using 900m but I switched it > in my example to my actual value (since I am going to use 300m for this > initial attempt at this). I will need to make it a lot more complicated > eventually! > > > > > On Sun, Nov 27, 2011 at 11:41 PM, Jim Lemon wrote: > >> On 11/28/2011 04:06 PM, Katrina Bennett wrote: >> >>> Sorry for not being more clear. >>> >>> I'll try to explain again. >>> >>> I have a rather large DEM and I need to scale daily temperature values >>> for >>> 10 years. >>> >>> I am using the sapply function to adjust temperatures (t.mean.1.c) based >>> on >>> lapse rates across DEM points (cdem) which meet the condition of >>> elevation >>> being greater than 300m. Below ~300m the lapse rate calculate changes >>> because temperatures decrease with elevation. Above ~300m the lapse rates >>> calculation must account for an inversion, where temperature increases >>> with >>> elevation. >>> >>> What I would like as a result is the values of the DEM which are lower >>> than >>> 300m to retain a generic lapse rate adjustment, and the DEM values great >>> than 300 to be treating using a different function. >>> >>> ... >>> >> Hi Katrina, >> If cdem is a matrix (and it looks like a vector in your example), you >> could try: >> >> le300<-function(elevation,**temp) return(temp+elevation*0.0065) >> gt300<-function(elevation,**temp) return(temp-elevation*0.0065) >> cdem[cdem<=300]<-le300(cdem[**cdem<=300],t.mean.1.c[cdem<=**300]) >> cdem[cdem>300]<-gt300(cdem[**cdem>300],t.mean.1.c[cdem>300]**) >> >> where le300 is the correction for altitudes less than or equal to 300M >> and gt300 is the one for altitudes greater than 300M. This works for me >> with a toy function. However, I don't get the same numbers as you and >> wonder if your function is doing the same as the ones I use above. Also, >> why 300m in the text and 900m in the functions? Are we mixing up feet and >> meters? >> >> Jim >> > > > > -- > Katrina E. Bennett > PhD Student > University of Alaska Fairbanks > International Arctic Research Center > 930 Koyukuk Drive, PO Box 757340 > Fairbanks, Alaska 99775-7340 > 907-474-1939 office > 907-385-7657 cell > kebenn...@alaska.edu > > > Personal Address: > UAF, PO Box 752525 > Fairbanks, Alaska 99775-2525 > bennett.katr...@gmail.com > > -- Katrina E. Bennett PhD Student University of Alaska Fairbanks International Arctic Research Center 930 Koyukuk Drive, PO Box 757340 Fairbanks, Alaska 99775-7340 907-474-1939 office 907-385-7657 cell kebenn...@alaska.edu Personal Address: UAF, PO Box 752525 Fairbanks, Alaska 99775-2525 bennett.katr...@gmail.com [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] Sum matrix by rows, conditional on value
I'd like to sum a matrix only where the matrix meets a specific condition. The matrix has 365 rows and about 50,000 columns. str(cdem.mat.yr) num [1:365, 1:41772] -43.5 -48.4 -45.9 -38.4 -32 ... I'm having trouble replicating this because it is not working out for me, so I'm unsure I can provide an solid working example (apologies). I would like to subset my matrix where the values are greater than zero. I can do this easily using the following command to generate a matrix of TRUE/FALSE conditions. thaw.index <- subset(cdem.mat.yr > 0) However, every time I then try to run apply, or rowSums to sum the matrix rows using this index, I get back errors messages or the wrong answer. There is a lot of values that meet this condition, so I know this issue is with my formatting of the argument. i.e. rowSums(cdem.mat.yr == thaw.index) [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [66] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [131] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [196] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [261] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [326] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 thaw.index.sum <- apply(cdem.mat.yr[which(cdem.mat.yr > 0)], 1, sum) Error in apply(cdem.mat.yr[which(cdem.mat.yr > 0)], 1, sum) : dim(X) must have a positive length This just provides me with one single value (I want the sum of all the rows individually) sum(apply(cdem.mat.yr, 1, function(x) x>0)) If I only run the apply over the entire matrix, this works but I can't subset it according to my condition. apply(cdem.mat.yr, 1, sum) [1] -1834000.521 -2038684.652 -1932676.585 -1619369.113 -1353598.426 -1190377.640 -1263541.796 -1438280.178 -1472440.385 -1714465.774 [11] -1945920.377 -2163888.712 -1836656.208 -1772994.790 -1864650.604 -1633144.043 -1580619.187 -1684046.620 -1769963.843 -1787676.116 [21] -1643345.342 -1497455.795 -1580307.433 -1483559.628 -1531067.546 -1557093.271 -1363626.528 -1160882.203 -1118893.494 -1352922.958 [31] -1441715.250 -1539084.024 -1717835.433 -1806727.136 -1887120.912 -1645721.673 -1310700.520 -1531408.644 -1582011.715 -1460538.996 [41] -1192880.165 -121.092 -1190390.732 -1241594.368 -1180704.394 -1025346.594 -1073581.060 -889662.991 -798735.050 -784632.707 [51] -874676.120 -957186.890 -1054610.980 -1067208.121 -1222859.353 -1477327.599 -1635653.310 -1696308.770 -1473630.951 -1283105.867 [61] -1061390.704 -811017.224 -875804.422 -851608.811 -948160.325 -1440351.359 -1206523.958 -1143659.246 -1405071.144 -1421438.254 [71] -1374929.105 -1336184.952 -1237185.588 -1082307.120 -1019742.616 -958257.706 -888078.311 -790481.841 -821010.686 -907205.025 [81] -966761.676 -926937.928 -908394.310 -976085.444 -971934.490 -703952.655 -521541.649 -625973.624 -743458.875 -631452.421 [91] -584709.631 -565843.210 -604252.152 -616885.977 -522011.655 -576824.263 -726170.003 -822902.735 -897385.940 -668897.194 [101] -525227.323 -493291.723 -559480.809 -627790.133 -607923.974 -535240.664 -346627.878 -343257.607 -287171.179 -324723.615 [111] -389052.208 -420393.385 -498589.819 -542040.688 -394442.745 -183163.637 -126540.029 -186213.012 -179799.971 -364410.639 [121] -309555.880 -357052.251 -321362.137 -394878.460 -498785.071 -309942.686 -276417.534 -337700.381 -304804.510 -238100.600 [131] -261210.843 -201821.616 -299377.673 -232015.614 -121752.676 -154925.661 -145809.72915840.738 145755.754 -33601.212 [141] -24323.63035036.73155156.44148603.82487203.646 139653.449 111722.558 101036.307 153884.464 153151.263 [151] 112680.914 108730.812 110198.055 127087.03377174.238 -67632.638 -35129.97656801.006 6712.631 8838.200 [161]40086.874 -29691.225 -55861.564 7561.50491232.944 31752.447 1694.756 -43835.544 2522.88341727.218 [171]26918.990 128692.011 114752.327 131455.862 127149.113 144686.214 160344.465 102204.088 130322.78570392.818 [181] 100384.523 210138.826 230235.443 211137.595 199770.103 167185.988 103951.789 125589.859 224791.286 302672.172 [191] 268148.251 258709.018 263356.469 122460.767 103103.124 1996.53055150.667 148763.608 188425.704 172693.650 [201] 173253.653 101070.947 142112.846 112317.201 101550.125 157215.618 184009.18360265.75094310.493 123499.949 [211] 174061.906 247635.961 268679.388 277766.454 307758.389 357310.030 178755.543 127887.6049997
Re: [R] Sum matrix by rows, conditional on value
Yes, that worked for me. Thank you. On Mon, Nov 28, 2011 at 10:16 PM, David Winsemius wrote: > > On Nov 29, 2011, at 1:08 AM, Katrina Bennett wrote: > > I'd like to sum a matrix only where the matrix meets a specific condition. >> The matrix has 365 rows and about 50,000 columns. >> > > If you describe the meaning attached to the data if might help readers > understand what the acceptable options might be. (Please read the Posting > Guide.) > > > >> str(cdem.mat.yr) >> num [1:365, 1:41772] -43.5 -48.4 -45.9 -38.4 -32 ... >> >> I'm having trouble replicating this because it is not working out for me, >> so I'm unsure I can provide an solid working example (apologies). >> >> I would like to subset my matrix where the values are greater than zero. I >> can do this easily using the following command to generate a matrix of >> TRUE/FALSE conditions. >> >> thaw.index <- subset(cdem.mat.yr > 0) >> > > Are you hoping to subset particular rows? > > If not, then you can make a copy and then set the negative values to NA > > thaw.mat <-cdem.mat.yr > is.na(thaw.mat)<- thaw.mat < 0 > > and ... > > rowSums(thaw.mat, na.rm=TRUE) > > > >> However, every time I then try to run apply, or rowSums to sum the matrix >> rows using this index, I get back errors messages or the wrong answer. >> > > I am guessing that your subset operation produced a vector and that apply > or rowSums no longer makes any sense. > > There is a lot of values that meet this condition, so I know this issue is >> with my formatting of the argument. >> >> i.e. >> rowSums(cdem.mat.yr == thaw.index) >> [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> [66] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> [131] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> 0 >> 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> [196] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> 0 >> 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> [261] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> 0 >> 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> [326] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> 0 >> 0 0 0 0 0 >> >> thaw.index.sum <- apply(cdem.mat.yr[which(cdem.**mat.yr > 0)], 1, sum) >> Error in apply(cdem.mat.yr[which(cdem.**mat.yr > 0)], 1, sum) : >> dim(X) must have a positive length >> >> This just provides me with one single value (I want the sum of all the >> rows >> individually) >> sum(apply(cdem.mat.yr, 1, function(x) x>0)) >> >> If I only run the apply over the entire matrix, this works but I can't >> subset it according to my condition. >> >> apply(cdem.mat.yr, 1, sum) >> [1] -1834000.521 -2038684.652 -1932676.585 -1619369.113 -1353598.426 >> -1190377.640 -1263541.796 -1438280.178 -1472440.385 -1714465.774 >> [11] -1945920.377 -2163888.712 -1836656.208 -1772994.790 -1864650.604 >> -1633144.043 -1580619.187 -1684046.620 -1769963.843 -1787676.116 >> [21] -1643345.342 -1497455.795 -1580307.433 -1483559.628 -1531067.546 >> -1557093.271 -1363626.528 -1160882.203 -1118893.494 -1352922.958 >> [31] -1441715.250 -1539084.024 -1717835.433 -1806727.136 -1887120.912 >> -1645721.673 -1310700.520 -1531408.644 -1582011.715 -1460538.996 >> [41] -1192880.165 -121.092 -1190390.732 -1241594.368 -1180704.394 >> -1025346.594 -1073581.060 -889662.991 -798735.050 -784632.707 >> [51] -874676.120 -957186.890 -1054610.980 -1067208.121 -1222859.353 >> -1477327.599 -1635653.310 -1696308.770 -1473630.951 -1283105.867 >> [61] -1061390.704 -811017.224 -875804.422 -851608.811 -948160.325 >> -1440351.359 -1206523.958 -1143659.246 -1405071.144 -1421438.254 >> [71] -1374929.105 -1336184.952 -1237185.588 -1082307.120 -1019742.616 >> -958257.706 -888078.311 -790481.841 -821010.686 -907205.025 >> [81] -966761.676 -926937.928 -908394.310 -976085.444 -971934.490 >> -703952.655 -521541.649 -625973.624 -743458.875 -631452.421 >> [91] -584709.631 -565843.210 -604252.152 -616885.977 -522011.655 >> -576824.263 -726170.003 -822902.735 -897385.940 -668897.194 >> [101] -525227.323 -493291.723 -559480.809 -627790.133 -607923.974 >> -535240.664 -346627.878 -343257.607 -287171.179 -324723
[R] Find position of asymptote
Hi all,
I would like to find the x position of an two asymptotes.
Here is a sample of what I would like to do:
x <- seq(1, 153,, 153)
a <- 85
m <- 65
s =-1.5
fn <- function (x, a, m, s) { a * (exp((m - x)/s) * (1/s))/((1 +
exp((m - x)/s)))^2 }
plot.deriv1 <- fn(1:153, a, m, s)
I can find the midpoint/minimum of this function easily. However, I
can not find the start and end of the two asymptote position (along
the x).
mid.point <- optimize(f=fn.sslogis.deriv1, interval=c(1:153), a=a,
m=m, s=s)$minimum
Thanks for your assistance,
Katrina
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Find position of asymptote
Hi David, thanks for the reply.
This is not a homework problem, although it may sound like one :) I
was trying to provide a reproducible example of what I am trying to
do.
The problem is something I am trying to work on for my PhD program.
I've been using the nls() function to derive a self-starting logistic
function to calculate snowmelt duration in Alaskan watersheds.
Here's a snippet of the nls function calculation **not reproducible**
#year 2000
print("year 2000")
Asym.2000 <- mean(na.omit(dat.cl.2000$dat)[1:10])
na.dat.2000 <- which(!is.na(dat.cl.2000$dat))
loess.smooth.2000 <- loess.smooth(dat.cl.2000$x[na.dat.2000],
dat.cl.2000$dat[na.dat.2000], span=0.20)
diff.loess.smooth.2000 <- which(diff(loess.smooth.2000$y) ==
(min(diff(loess.smooth.2000$y))), arr.ind=T)
val.in.jday.2000 <- loess.smooth.2000$y[diff.loess.smooth.2000]
xmid.2000 <- val.in.jday.2000
scal.2000 <- (min(na.omit(dat.cl.2000$dat)) -
max(na.omit(dat.cl.2000$dat))) / (max(na.omit(dat.cl.2000$x)) -
min(na.omit(dat.cl.2000$x))) * 10
fit.dat.2000 <- nls(dat ~ SSlogis(x, Asym, xmid,scal), data =
dat.cl.2000, start=list(Asym=round(Asym.2000), xmid=round(xmid.2000),
scal=scal.2000), control = list(maxiter = 500, warnOnly = TRUE),
trace=TRUE)
What I want to be able to do however, is to determine the start and
end points (dates) of the asymptote of the curve when the snowmelt
initiates.
Prior to this I was using the second derivative mins and maxes.
#init <- optimize(f=fn.sslogis.deriv2, interval=c(1:mid.point), a=a,
m=m, s=s)$minimum #this is the start of the snowmelt season
#term <- optimize(f=fn.sslogis.deriv2, interval=c(mid.point:153), a=a,
m=m, s=s, maximum=TRUE)$maximum #this is the end of the snowmelt
season
However, this cuts the melt too short.
Any assistance on this non-hw problem would be appreciated.
Sincerely,
Katrina
On Thu, Apr 19, 2012 at 4:34 AM, David Winsemius wrote:
>
> On Apr 19, 2012, at 4:41 AM, Katrina Bennett wrote:
>
>> Hi all,
>>
>> I would like to find the x position of an two asymptotes.
>>
>> Here is a sample of what I would like to do:
>>
>> x <- seq(1, 153,, 153)
>> a <- 85
>> m <- 65
>> s =-1.5
>> fn <- function (x, a, m, s) { a * (exp((m - x)/s) * (1/s))/((1 +
>> exp((m - x)/s)))^2 }
>> plot.deriv1 <- fn(1:153, a, m, s)
>>
>> I can find the midpoint/minimum of this function easily. However, I
>> can not find the start and end of the two asymptote position (along
>> the x).
>
>
> This is not a homework help-line.
>
>
>> mid.point <- optimize(f=fn.sslogis.deriv1, interval=c(1:153), a=a,
>> m=m, s=s)$minimum
>>
>> Thanks for your assistance,
>>
>> Katrina
>
>
> --
>
> David Winsemius, MD
> West Hartford, CT
>
--
Katrina E. Bennett
PhD Student
University of Alaska Fairbanks
International Arctic Research Center
930 Koyukuk Drive, PO Box 757340
Fairbanks, Alaska 99775-7340
907-474-1939 office
907-385-7657 cell
kebenn...@alaska.edu
Personal Address:
UAF, PO Box 752525
Fairbanks, Alaska 99775-2525
bennett.katr...@gmail.com
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] Find position of asymptote
Dear R Help,
Sorry I wasn't more clear before. Here is another crack at this.
What I am still trying to do is estimate the point on a line when the
slope changes or asymptotes. I have found some similar postings
talking about this but no answers.
https://stat.ethz.ch/pipermail/r-help/2003-January/028847.html for example
http://stackoverflow.com/questions/8245792/r-marking-slope-changes-in-loess-curve-using-ggplot2
(best one I see)
Based on the above, I've modified my code. I have provided below some
sample data so you can see the entire chain of work that summarizes
the issue I am trying to resolve.
dat <-
c(75.44855206,NA,NA,NA,82.70745342,82.5335019,88.56617647,80.00128866,94.15418227,86.63987539,93.91052952,74.10612245,86.62289562,90.47961047,NA,NA,82.45320197,72.14371257,NA,71.44104803,72.59742896,68.36363636,NA,NA,61,NA,NA,71.26502909,NA,85.9333,84.34248284,79.00522193,79.64223058,97.2074017,88.43700548,96.40413877,95.13511869,92.57379057,93.97498475,NA,97.55995131,89.53321146,97.21728545,93.21980198,77.54054054,95.85392575,86.25684723,97.55325624,80.03950617,NA,91.34023128,92.42906574,88.59433962,65.77272727,89.63772455,NA,NA,NA,NA,74.86344239,83.57594937,70.22516556,65.30543319,NA,NA,67.84852294,60.90909091,54.79303797,NA,52.18735363,33.47003155,NA,41.34693878,24.5047043,NA,NA,NA,NA,9.9,13.6875,NA,11.90267176,84.14285714,3.781456954,NA,1.432926829,4.26557377,1.823529412,0.444620253,4.711155378,NA,6.320284698,0.581632653,0.144578313,3.7,0,0,0,0,0,NA,0.032947462,0,0,10.54545455,0,NA,0.561007958,0.75,NA,0.048780488,0.74137931,NA,2.023339318,0,0,0!
,NA,NA,0.156950673,NA,0.283769634,32.81818182,NA,NA,0,NA,0,0,0,NA,0.212454212,3.120181406,NA,0.011811024,NA,0,0.120430108,5.928571429,1.75,0.679292929,0.97,NA,0,NA,NA,1,0.38547486,NA,1.460732984,0.007795889,0.05465288,0.004341534)
x.seq <- seq(1, 153,, 153)
dat.cl.2000 <- as.data.frame(cbind(x.seq, dat))
names(dat.cl.2000) <- c("x", "dat")Asym.2000 <-
mean(na.omit(dat.cl.2000$dat)[1:10])
na.dat.2000 <- which(!is.na(dat.cl.2000$dat))
loess.smooth.2000 <- loess.smooth(dat.cl.2000$x[na.dat.2000],
dat.cl.2000$dat[na.dat.2000], span=0.20)
diff.loess.smooth.2000 <- which(diff(loess.smooth.2000$y) ==
(min(diff(loess.smooth.2000$y))), arr.ind=T)
val.in.jday.2000 <- loess.smooth.2000$y[diff.loess.smooth.2000]
xmid.2000 <- val.in.jday.2000
scal.2000 <- (min(na.omit(dat.cl.2000$dat)) -
max(na.omit(dat.cl.2000$dat))) / (max(na.omit(dat.cl.2000$x)) -
min(na.omit(dat.cl.2000$x))) * 10
fit.dat.2000 <- nls(dat ~ SSlogis(x, Asym, xmid,scal), data =
dat.cl.2000, start=list(Asym=round(Asym.2000), xmid=round(xmid.2000),
scal=scal.2000), control = list(maxiter = 500, warnOnly = TRUE),
trace=TRUE)
year.ind <- 2000
fit.dat <- get(paste("fit.dat.", year.ind, sep=""))
a <- coef(fit.dat)[[1]]
m <- coef(fit.dat)[[2]]
s <- coef(fit.dat)[[3]]
sslogis <- expression(a/((1+exp((m-x)/s
sslogis.deriv1 <- D(sslogis, "x")
sslogis.deriv2 <- D(D(sslogis, "x"), "x")
fn.sslogis <- function (x, a, m, s) { a/((1+exp((m-x)/s))) }
fn.sslogis.deriv1 <- function (x, a, m, s) { a * (exp((m - x)/s) *
(1/s))/((1 + exp((m - x)/s)))^2 }
fn.sslogis.deriv2 <- function (x, a, m, s) {
-(a * (exp((m - x)/s) * (1/s) * (1/s))/((1 + exp((m - x)/s)))^2 -
a * (exp((m - x)/s) * (1/s)) * (2 * (exp((m - x)/s) * (1/s) *
(1 + exp((m - x)/s)/1 + exp((m - x)/s)))^2)^2) }
asym <- a
slope <- s
mid.point <- optimize(f=fn.sslogis.deriv1, interval=c(1:153), a=a,
m=m, s=s)$minimum #this will find the minimum point halfway through
melt
init <- optimize(f=fn.sslogis.deriv2, interval=c(1:mid.point), a=a,
m=m, s=s)$minimum #this is the start of the snowmelt season
term <- optimize(f=fn.sslogis.deriv2, interval=c(mid.point:153), a=a,
m=m, s=s, maximum=TRUE)$maximum #this is the end of the snowmelt
season
duration <- init-term
~
Ok, this part of the code is where I really want to try and calculate
the points (init and term) where the curve asymptotes (using the first
derivative). On the below plot, this should be at ~40 and ~110 (these
represent days in my analysis).
plot.deriv1 <- fn.sslogis.deriv1(x.seq, a, m, s)
plot(plot.deriv1)
I found in the help pages above that I could do this based only on the
loess curve. I can take the diff of the loess prediction to calculate
a new function that has two local minimums.
loess.2000 <- loess(dat.cl.2000$x[na.dat.2000]~dat.cl.2000$dat[na.dat.2000])
loess.predict <- predict(loess.2000, newdata=x.seq)
loess.predict1 <- diff(loess.predict)
loess.predict2 <- c(NA,loess.predict1)
plot(x.seq, loess.predict2, main="diff(loess model)")
However, this still isn't providing me a means to get the points of
transition along this line.
Is there a better way to then pick off the change points or find the
asymptotes of a function in R?
Thank you.
Katrina
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do r
