[R] Extracting tolerance in R?
Dear list, How is the tolerance for a model parameter in an lm() call extracted? I did not see a solution in the documentation for lm(), or predict(), nor in the archives using 'tolerance' as the search string. I also checked into the nlme package, though nothing popped out at me. Sincerely, KeithC. __ R-help@stat.math.ethz.ch 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] {nuSpectral} irreg sampled ts, phase based periodogram
Dear List-mates,
I just read Mathias, et al (2004) about their {nuSpectral} package in R for
a weighted leased squares method for unevenly sampled time-series.
What kind of experiences have users had with the package?
R-Site Search came up blank. I've never asked about a package that was
downloaded outside of CRAN. Feel free to redirect me as the case may be.
Mathias, Grond, Guardans, Seese, Canela, et al (2004). Algorithms for
spectral analysis of irregularly sampled time series. Journal of Statistical
Software. 11-2.
Rgds,
KeithC.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] MenuRead() Question
Dear List-mates,
I think the difficulty I'm having is localized to the MenuType() call made
from within MenuRead(). I'm not used to seeing text operations such as ^
or [], so am having trouble understanding what's going on.
I'm interested in understanding the regular expression:
regexpr(^[$]Tk[.].+/, menu), and why it my menus.txt file is not returning
at that point.
?regex explains that the ^ symbol excludes the text enclosed in brackets
what is in the character string in following brackets, so a bunch of
matching on string vectors going on that I don't understand well enough yet.
Do I need to install PCRE-6.6 for this to start working?
Rgds,
KeithC.
-Original Message-
From: Keith Chamberlain [mailto:[EMAIL PROTECTED]
Sent: Saturday, March 25, 2006 4:26 PM
To: 'r-help@stat.math.ethz.ch'
Cc: 'Philippe Grosjean'
Subject: MenuRead() Question
Dear List-mates,
I'm trying to read a tk window menu from a file using {svWidgets} and
'menus.txt' but am receiving Warnings without seeing the desired
consequences of the call.
library(svWidgets)
tkWinAdd(KSesnMain,title=kLab Session Manager for R, pos=+0+0)
MenuRead(file=menus.txt)
Warning messages:
1: Unrecognized menu type for $KSesnMain/MenuTest in: MenuType(menu)
2: Unrecognized menu type for $KSesnMain/MenuTest in: MenuType(menu)
3: Unrecognized menu type for $KSesnMain/MenuTest in: MenuType(menu)
4: Unrecognized menu type for $KSesnMain/MenuTest in: MenuType(menu)
menu.txt - in RHOME
$KSesnMain
|$MenuTest
||Objects ~~ ls()
||-
||Path ~~ search()
Please Advise,
KeithC.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] MenuRead() Question
Dear Gabor Philipe,
Thank you for clarifying regular expressions, and how tk menus are
demarcated in menus.txt. Regular expressions are beginning to make sense for
me, and the menu is now being added properly to the tk window.
regexpr(^[$]Tk[.].+/, menu)...
If I understand better, the expression is matching $Tk. but since $ and
. have special meanings in regular expressions, they need to be in
brackets so that they are interpreted literally.
Rgds,
KeithC.
-Original Message-
From: Philippe Grosjean [mailto:[EMAIL PROTECTED]
Sent: Sunday, March 26, 2006 9:49 AM
To: Keith Chamberlain
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] MenuRead() Question
Hi Keith,
If you want to define a Tk menu using MenuRead() and a text file to
define your menu, you have to start it with Tk.. That way, MenuRead()
recognizes that it is a Tk menu. So, rewrite your menu definition file as:
menu.txt - in RHOME
$Tk.KSesnMain
|$MenuTest
||Objects~~ ls()
||-
||Path ~~ search()
Note also that KSesnMain must point to a valid Tk window previously
constructed.
Best,
Philippe Grosjean
Keith Chamberlain wrote:
Dear List-mates,
I think the difficulty I'm having is localized to the MenuType() call made
from within MenuRead(). I'm not used to seeing text operations such as ^
or [], so am having trouble understanding what's going on.
I'm interested in understanding the regular expression:
regexpr(^[$]Tk[.].+/, menu), and why it my menus.txt file is not
returning
at that point.
?regex explains that the ^ symbol excludes the text enclosed in brackets
what is in the character string in following brackets, so a bunch of
matching on string vectors going on that I don't understand well enough
yet.
Do I need to install PCRE-6.6 for this to start working?
Rgds,
KeithC.
-Original Message-
From: Keith Chamberlain [mailto:[EMAIL PROTECTED]
Sent: Saturday, March 25, 2006 4:26 PM
To: 'r-help@stat.math.ethz.ch'
Cc: 'Philippe Grosjean'
Subject: MenuRead() Question
Dear List-mates,
I'm trying to read a tk window menu from a file using {svWidgets} and
'menus.txt' but am receiving Warnings without seeing the desired
consequences of the call.
library(svWidgets)
tkWinAdd(KSesnMain,title=kLab Session Manager for R, pos=+0+0)
MenuRead(file=menus.txt)
Warning messages:
1: Unrecognized menu type for $KSesnMain/MenuTest in: MenuType(menu)
2: Unrecognized menu type for $KSesnMain/MenuTest in: MenuType(menu)
3: Unrecognized menu type for $KSesnMain/MenuTest in: MenuType(menu)
4: Unrecognized menu type for $KSesnMain/MenuTest in: MenuType(menu)
menu.txt - in RHOME
$KSesnMain
|$MenuTest
||Objects ~~ ls()
||-
||Path ~~ search()
Please Advise,
KeithC.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] MenuRead() Question
Dear List-mates,
I'm trying to read a tk window menu from a file using {svWidgets} and
'menus.txt' but am receiving Warnings without seeing the desired
consequences of the call.
library(svWidgets)
tkWinAdd(KSesnMain,title=kLab Session Manager for R, pos=+0+0)
MenuRead(file=menus.txt)
Warning messages:
1: Unrecognized menu type for $KSesnMain/MenuTest in: MenuType(menu)
2: Unrecognized menu type for $KSesnMain/MenuTest in: MenuType(menu)
3: Unrecognized menu type for $KSesnMain/MenuTest in: MenuType(menu)
4: Unrecognized menu type for $KSesnMain/MenuTest in: MenuType(menu)
menu.txt - in RHOME
$KSesnMain
|$MenuTest
||Objects ~~ ls()
||-
||Path ~~ search()
Please Advise,
KeithC.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Platform independent dialogs menus?
Dear list mates,
Are {utils} dialog box functions, and winMenuAdd... functions used to change
(e.g. Console) menus, platform dependent?
I'm writing a script loaded with .First that provides first time users in a
lab course with the ability to select, load, and change between what I
called a 'session' (more specific save that focuses on the session object I
defined rather than the whole workspace, intended to run many different
sessions through the course of what would be one workspace).
I'm using winMenuAdd() calls to generate their 'Session' menu at startup,
the menus call functions sourced in for the menu actions. In the menu
actions, I call routines that use select.list() and file.choose() calls to
interact with users.
I do not work with Macs often, and from what I've gathered today in posts
about cross-platform difficulties, my sense of being intimidated seems
well placed to me (then again, breaks some sleep would probably help). I
have not had the chance to test routines on a Mac yet, so I have no idea
what to expect. Is this tract I took with winMenuAdd() related [{utils}
windows build] an appropriate route to take wrt the Mac build of R, or would
I be better off using another package?
Please advise,
KeithC.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Platform independent dialogs menus?
Dear Philippe, ( list)
I am not sure what is meant by 'floating window' but thank you for
clarifying the 'win' prefix. I had been considering different possible
senses of that prefix: win as in 'ts' window windowing functions (clearly
not), a generic 'win'dow function (abbreviated) for some GUI, or ms windows
GUI windowing function in particular. I'm noticing that the description in
{utils:winMenuAdd} indicates 'for windows', and it's quite clear (now,
anyway) that it is to be taken as ms windows. But cool, I WAS barking down
the wrong tree can now refocus efforts. Curious how the correct meaning
did not pop out at me.
I was not aware of the {Session} package. What I've merely glanced at
suggests the package is exactly (or really, really darn close) to what I had
in mind. I appreciate the suggestion.
I downloaded SciViews Tinn-R 2 days ago, and LOVE Tinn-R so far; syntax
highlighting, oh yes! My rudimentary understanding, however, was that
SciViews was only compiled for MS Windows so far. Can I run SciViews
natively on a Mac now, or is the cross-platform part specific to some of the
routines?
Respecftully,
KeithC.
-Original Message-
From: Philippe Grosjean [mailto:[EMAIL PROTECTED]
Sent: Monday, March 20, 2006 6:18 AM
To: Keith Chamberlain
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] Platform independent dialogs menus?
Keith Chamberlain wrote:
Dear list mates,
Are {utils} dialog box functions, and winMenuAdd... functions used to
change (e.g. Console) menus, platform dependent?
Yes, Windows-only as the 'win' prefix suggests it.
I'm writing a script loaded with .First that provides first time users in
a lab course with the ability to select, load, and change between what I
called a 'session' (more specific save that focuses on the session object
I defined rather than the whole workspace, intended to run many
different sessions through the course of what would be one workspace).
See also the 'session' package on CRAN for that.
I'm using winMenuAdd() calls to generate their 'Session' menu at startup,
the menus call functions sourced in for the menu actions. In the menu
actions, I call routines that use select.list() and file.choose() calls to
interact with users.
I do not work with Macs often, and from what I've gathered today in posts
about cross-platform difficulties, my sense of being intimidated seems
well placed to me (then again, breaks some sleep would probably help). I
have not had the chance to test routines on a Mac yet, so I have no idea
what to expect. Is this tract I took with winMenuAdd() related [{utils}
windows build] an appropriate route to take wrt the Mac build of R, or
would
I be better off using another package?
For a platform-independent way of defining menus (you will have a
floating window with your menu), look at ?MenuAdd in package svWidgets
(SciViews bundle). With these functions, you even have more control on
the menus (define shortcuts, trigger the menus through R code, for
instance), and you can define your menu in one R instruction and a
simple text file to describe the menu structure, like this:
# Create a Tk window then add a menu to it
$MyTkWindow
|$MenuTest
||Objects ~~ ls()
||-
||Path ~~ search()
# Add menus to the RGui console (Windows only)
$ConsoleMain
|$Testit
||One ~~ cat(One triggered!\n)
||-
||Two ~~ cat(Two triggered!\n) ~~ state = disable
||$SubMenu
|||Three ~~ cat(Three triggered!\n)
|||Four ~~ cat(Four triggered!\n)
||Five~~ cat(Five triggered!\n)
# Add menu to the RGui console popup (Windows only)
$ConsolePopup
|$TestitPopup
||Six ~~ cat(Six triggered!\n)
If the preceeding menu definition is in a file named Menus.txt in the
/gui subdirectory of your MyPackage package, you can do:
library(svWidgets)
tkWinAdd(MyTkWindow, title = Menu window, pos =-40+20)
MenuReadPackage(MyPackage)
... and you got all your menus configured at once (for MyTkWindows +
RGui console + RGui console popup menu
Best,
Philippe Grosjean
Please advise,
KeithC.
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] GAM using R tutorials?
Dear Templ ( Others), Thank you for pointing out the contributed papers. It is one of those things that I have been introduced to more than once, but I personally never remembered to check the section out while in the acute phase of problem-can't do it-what the heck is going on. Dear Michael, Regarding frustrating replies, those will happen. I lost sleep over such things in the past, only to learn that I didn't really have time to get caught up. It is normal when investing time and effort into a reply to want some kind of assurance that some legwork has been done. I think Gavin just didn't get the sought after assurance from your message. The reply may not have been particularly useful to you, per say, but it did serve as a reminder to a wider audience of subscribers who read posts. Given that this is done over email, and given the huge diversity of subscribers around the globe with different language histories, the error rate in communication could only increase from that of your or my locale. Hence, the importance of the posting guide FAQs. More often than not, however, I've later found replies that I previously considered terse or inappropriate to be quite accurate. On the other side, I walked through an RD lab once pointed out a solution to a router configuration problem that the technician had been working on for a few days. Having been trained on that one system, it was quite obvious what needed to be fixed. I expected the technician to be grateful, but he never spoke to me again! Well, that's my 2 cents. I do acknowledge that it wasn't requested. My hope is that someone out there in listland would find it useful, even if not you. Yes, I guess I am having a conversation with the wind... Rgds, KeithC. Thu, 16 Mar 2006 00:26:46 -0800 Michael [EMAIL PROTECTED] Re: [R] GAM using R tutorials? Hi Templ, That's very helpful! Indeed I've printed out the gam portion and I am digesting now... Thank you so much! I really appreicate your constructive and helpful advice! Best, Michael. On 3/15/06 TEMPL Matthias [EMAIL PROTECTED] wrote: Have you looked at: An Introduction to R: Software for StatisticalModelling Computing by Petra Kuhnert and Bill Venables which is available at http://cran.r-project.org/other-docs.html Hope this helps. Best, Matthias Thu, 16 Mar 2006 00:25:56 -0800 Michael [EMAIL PROTECTED] Re: [R] GAM using R tutorials? Interesting! Very interesting! You seem to assume that I haven't done anything as you've listed below: I have even listed a book that I was aware of but I could not find in our local library... and doing some research and I have already started playing with gam, and everybody on earth knows that you can do ?gam and help.search() in R to find info, but after I have done all those, I still feel that I need some more infohelp... I don't know how you reach your ungrounded conclusion? I think you are wasting everybody's bandwidth and time by not providing constructive suggestions but discouraging (and overly repetitive) comments disregarding the effort a newbie had already paid... I think you really need to consult posting guide yourself... And remember, you can always remain silent -- you don't need to show off that you are an experienced user and for newbie your first response is always something like scolding: why didn't you search? I strongly believe that your answer is inappropriate!!! __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Trellis - setting xlim or ylim by data range in whole column or row
Dear List-mates, I have been trying to set up a 4x8 trellis plot (that is, 4 columns, 8 rows), and would like to have the axis limits set based on the data range of rows (for ylim) and columns (for xlim). I've been using the call: foo-xyplot(y~x|Epoch+Subject, type=c(l,r), par.strip.text=list(cex=0.5), ...) and updating to see different effects of scale adjustments, etc. foo-update(foo, scale=list(relation=sliced)) #etc. I can have each panel with its own limits using relation=free, limits wide enough to accommodate the whole data range with same or the limits of the widest panel for all panels with split. I have not, however, figured out a way to have the limits for x y grouped by their column or row. Documentation points to 3 alternate ways of trying to set stuff in panels. (1) a prepanel function, (2) using scales, or (3) focusing in on each panel individually changing some setting. I've played around with accessing different panels individually with foo[column, row], and using a list to determine which get displayed (instead of using skip because I can't get skip to work). Would I be able to set the xlim values in a similar way? foo[1,]$ylim-newvalues to set a whole columns ylims (e.g by data range of y in conditioning variable 2 (subject in my case)) and foo[,1]$xlim-newvalues to get a whole rows xlims by the data range of x in each level of conditioning variable 1 (Epoch in my formula))? If so, what attribute should I access, and if not what would you recommend? I've been reading posts, working examples in Venables Ripley 4th Ed., and experimenting with different things for the last 4 days. I'm still not used to the lattice terminology, so I could have easily miss interpreted what something was meant for (example- prepanel makes no sense to me yet). Conversely, I got a lot farther, a lot faster, using Lattice than I did using plot or ts.plot. In addition to a much shorter list of attributes that don't make sense to me yet than otherwise, I have been really tickled with the Lattice package. Thanks in advance for your feedback. KeithC. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] Trellis - setting xlim or ylim by data range in whole column or row
Dear Deepayan, My deepest thanks! Your example code worked perfectly. I understood the caveats you detailed about figuring out the axis limits before resolving the layout, which depend on a bunch of other stuff. I'm glad this is possible. Yes, this particular layout does only make sense in this special case (e.g. 2 conditioning variables exactly n and m-levels). Thank you again for your feedback. I am excited to go put your examples to work. Rgds, KeithC. -Original Message- From: Deepayan Sarkar [mailto:[EMAIL PROTECTED] Sent: Thursday, March 09, 2006 9:55 PM To: Keith Chamberlain Cc: r-help@stat.math.ethz.ch Subject: Re: Trellis - setting xlim or ylim by data range in whole column or row On 3/9/06, Keith Chamberlain [EMAIL PROTECTED] wrote: Dear List-mates, I have been trying to set up a 4x8 trellis plot (that is, 4 columns, 8 rows), and would like to have the axis limits set based on the data range of rows (for ylim) and columns (for xlim). I've been using the call: foo-xyplot(y~x|Epoch+Subject, type=c(l,r), par.strip.text=list(cex=0.5), ...) and updating to see different effects of scale adjustments, etc. foo-update(foo, scale=list(relation=sliced)) #etc. I can have each panel with its own limits using relation=free, limits wide enough to accommodate the whole data range with same or the limits of the widest panel for all panels with split. I have not, however, figured out a way to have the limits for x y grouped by their column or row. Documentation points to 3 alternate ways of trying to set stuff in panels. (1) a prepanel function, (2) using scales, or (3) focusing in on each panel individually changing some setting. I've played around with accessing different panels individually with foo[column, row], and using a list to determine which get displayed (instead of using skip because I can't get skip to work). Would I be able to set the xlim values in a similar way? foo[1,]$ylim-newvalues to set a whole columns ylims (e.g by data range of y in conditioning variable 2 (subject in my case)) and foo[,1]$xlim-newvalues to get a whole rows xlims by the data range of x in each level of conditioning variable 1 (Epoch in my formula))? If so, what attribute should I access, and if not what would you recommend? What you want is not impossible, but not simple either. There's a good reason for that---this behavior makes sense only in a very special situation: when you have exactly two conditioning variables (a and b, say) and your layout has exactly nlevels(a) columns and nlevels(b) rows. To fix ideas, consider the following example: library(lattice) d - expand.grid(a = gl(4, 1, 20), b = gl(8, 1, 40)) d$x - with(d, rnorm(nrow(d), mean = as.numeric(a) + as.numeric(b))) d$y - with(d, rnorm(nrow(d), mean = as.numeric(a) + as.numeric(b))) foo is as in your example, foo - xyplot(y ~ x | a + b, d) At this point, foo has no layout defined (foo$layout is NULL). The layout will be determined when it is plotted. Since there are exactly two conditioning variables, the default here is layout = dim(foo) = c(4, 8). So plotting foo is the same as update(foo, layout = dim(foo)) However, if you instead do update(foo, layout = c(0, prod(dim(foo (which is essentially what happens when there is one conditioning variable) you will get a different layout, probably 6x6 (unless you have resized the device). In this case, choosing limits by row or column no longer makes sense. You might say, why not do this whatever the layout, whether it makes sense or not. The problem is, the axis limits for each panel needs to be determined _before_ the layout, because the layout may depend on the aspect ratio and the aspect ratio may depend on the axis limits (e.g. when aspect = xy). Workarounds: It is possible to explicitly specify a limit for each panel as a list. For example, you could have (update doesn't work well for this): xyplot(y ~ x | a + b, d, scales = list(x = list(relation = free, limits = rep(list(c(0, 11), c(1, 12), c(2, 13), c(4, 14)), 8))), par.strip.text = list(cex = 0.5)) If you are happy with this, that's great, but you will most likely want to omit all but one set of labels and use the freed up space. To control the labelling, you can specify the tick mark locations as a list just like the limits. An undocumented trick to simplify this is to have TRUE for the defaults, and NULL to omit them, so you could do: xyplot(y ~ x | a + b, d, scales = list(x = list(relation = free, limits = rep(list(c(0, 11), c(1, 12), c(2, 13), c(4, 14)), 8), at = rep(list(TRUE, NULL), c(4, 28, par.strip.text = list(cex = 0.5)) Finally, to make use of the space: xyplot(y ~ x | a + b, d, scales = list(x = list(relation
[R] LME Correlation Component using spectrum()?
Dear List-mates, Is there a corStruct constructor already written to calculate the shape of the spectral density for linear models using Fourier estimates (e.g. terms for a linear model derived from the frequency domain)? I have data with a long memory process but do not want to destroy it by taking n-differences for a suitable corAR, or corARIMA object. Thank you, KeithC. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Re: [R] How do I normalize a PSD?
Dear Tom, Short answer, if your using spec.pgram(), use the smoothing kernel to get a better estimate at the frequency centered in the bandwidth. If your frequency bin of interest is wider than the bandwidth of the kernel, average across frequencies (I think). The estimate appears to be normalized already. If you are calculating your PSD independently, then oversample (e.g. 2, perhaps 4 or more times the highest frequency you need to characterize, and sum (not average) across the frequency bin to get the power estimate of the center frequency in that bin. Shove the whole signal into memory and don't window it if you can. The oversampling improves your estimate and reduces the variance. Getting the FFT of the whole signal at once eliminates edge effects that would occur if you segmented the data. Whatever your PSD is normalized to, the sum across frequencies in your bin will maintain that property. If you are really wanting to see the PSD plotted by period rather than frequency, plot the inverse of your frequency values. The normalization part is a separate (and huge, perhaps neglected, perhaps overrated) subject. Here is how I have tried to make sense out of it. I do not have years and years of experience in spectral analysis, but perhaps you may find this a more simple (minded?) description. I was taught the definition of the power spectrum using Numerical Recipes from Press, et al (e.g. any of 1988, 1992, and 2002, for FORTRAN, C, and C++ respectively). Their description seems different (as Leif suggested) from other sources. They also happen to caution (strongly) about the distinction between one-sided and two-sided PSDs, but I'm not sure to what extent other authors have concerned themselves with that point, or that it is even necessary in most practical applications (warning: non-expert opinion). In the definition I learned, a periodogram is considered normalized in the sense that the sum across the frequency range has something to do with some quantity (Sum Squared, Mean Squared, etc.) of a signal in the time domain. Let N = Number of samples Fc = Nyquist critical frequency C = Normalization parameter P2(f) = C * Mod(fft(signal))^2 from -fc to fc: Press et al call a two-sided estimate, even if just the power at positive frequencies is returned. P1(f) = C * Mod(fft(signal))^2 from DC to fc, but all frequencies excluding DC and Nyquist are multiplied by 2. Press et al. consider this the one-sided estimate. The reason they call those one or two sided is that the sum of the power over all of the frequencies will be equivalent even though the one sided estimate has only half the data points. Sum(P2) == sum(P1) The parameter C is what I think I understand as the normalization parameter, and besides the one/two sided nomenclature issue, is where authors seem to differ. Regardless of which definition of the PSD is used, where C = 1/N, the periodogram (to Press, et al.) would be normed to the Sum Squared of the signal in the time domain. I believe some authors might call that form un-normalized, depending on whether or not 1/N was included in the definition of a periodogram that they happen to have been working from. Where C = 1/N^2, the periodogram would be normed to the Mean Squared of the signal in the time-domain. Where C = some other factor, the periodogram is normed to the equivalent of that factor from the time-domain. That's my (perhaps oversimplified) understanding. From Venables Ripley, listed as a source on ?spec.pgram, the definition used for that function seems comparable - though I wouldn't count out the possibility that I missed or misread something. The factor C in that case is 1/N so the PSD estimate should sum to the Sum Squared in the time domain. There are 2 caveats. The first is a point that Leif pointed out, that the periodogram from spec.pgram() is normalized to frequency. I need to reread that section in Venables Ripley and look at the link provided in another recent reply in more detail to be certain. For now, I am making the assumption that the parameter would be used 'in addition' to the 1/N factor. I highly recommend that you reference the sources from the spec.pgram help page. Second, smoothing and tapering are defined separately from the PSD, so I do not think the equalities hold up as well between a quantity of the (unsmoothed) time-domain signal, and its normalized (smoothed /or tapered) power spectrum. I hope this helps, wasn't to long, and that my weak editing skills did not make to grand an appearance. Regards, KeithC. Message: 30 Date: Tue, 31 Jan 2006 13:02:03 -0800 From: Leif Kirschenbaum [EMAIL PROTECTED] Subject: Re: [R] How do I normalise a power spectral density To: r-help@stat.math.ethz.ch, Tom C Cameron [EMAIL PROTECTED] Message-ID: [EMAIL PROTECTED] Content-Type: text/plain; charset=iso-8859-1 I have done a fair bit of spectral analysis, and hadn't finished collecting my thoughts for a reply, so hadn't replied yet. What exactly do you mean
[R] Relating Spectral Density to Chi-Square distribution
Dear list, I had some confusion regarding what function too use in order too relate results from spec.pgram() too a chi-square distribution. The documentation indicates that the PSD estimate can be approximated by a chi-square distribution with 2 degrees of freedom, but I am having trouble figuring out how to do it in R, and figuring out what specifically that statement in the documentation means. I have very little exposure to distribution functions in R. I have started with a signal that is simply a sine function, with a visible periodicity within the nyquist range. a- sin(2*pi *(0:127)/16) plot(a, type=l) I then get the raw power spectrum using, for instance: PSD - spec.pgram(a, spans=NULL, detrend=F) Next I did a search for information on chi-square. help.search(chi square) . provided a list of potentials, of which Chisquare() seemed to fit because it mentioned the distribution. ?Chisquare . provided the documentation, which shows four functions and their descriptions. I've assumed that I need too use dchisq() for my purposes, so that the fitted distribution would be: [assuming the power returned by spec.pgram() ARE regarded as quantiles.] plot(dchisq(PSD$spec, df=2)) . but the values are not between (0,1). plot(PSD$freq, pchisq(PSD$spec, df=2, lower.tail=F)) . looks better because values range from 0-1. Please clarify how to fit the PSD estimate to a chi-square distribution and any transforms that may be needed to get the quantiles from the PSD. Is what I tried what the documentation 'had in mind' when it says df: The distribution of the spectral density estimate can be approximated by a chi square distribution with 'df' degrees of freedom. ? If so, what is the appropriate function call too use (pchisq(), dchisq(), or qchisq())? If not, what function should I consider? Thanks in advance, Keith C. [[alternative HTML version deleted]] __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] spec.pgram() normalized too what?
Dear list, What on earth is spec.pgram() normalized too? If you would like to skip my proof as to why it's not normed too the mean squared or sum squared amplitude of the discrete function a[], feel free too skip the rest of the message. If it is, but you know why it's not exact in spec.pgram() when it should be, skip the rest of this message. The issue I refer herein refers only too a single univariate time series, and may not reflect the issues encountered in 3 or more dimensions. I've been using Numerical Recipes too confirm my own routines, because the periodogram estimate is not normalized to the sum squared, nor the mean squared of the original discrete function, a[]. I'd expect the non-overlapped, and non-tapered version of the periodogram too sum to EXACTLY too some normalization of the discrete signal a[], and the word estimate should come into play only when making inferences about the real signal 'a' that the discretely sampled signal came from. a - sin(2*pi*(0:127)/16) N - length(a) # 128 PSD- spec.pgram(a, spans=NULL, detrend=F) ## Sum Squared amplitude of a[t] on [0, 127]. sum(abs(a)^2) [1] 64 ## Mean Squared amplitude of a[t] on [0, 127] sum(a^2)/N [1] 0.5 By Parseval's theorem, the integral of the one-sided PSD over zero and positive frequencies should equal the mean squared or sum squared amplitude of the discrete signal a[], assuming the PSD is normalized too the mean-square or sum squared of the signal a[] respectively. ## Integral of the PSD returned by spec.pgram() does not equal MS or SS of a[] sum(PSD$spec) [1] 32.28501 Since the integral over positive frequencies does not equal the mean or sum squared of the signal, I did some digging. The documentation for spec.pgram() states that only the power at positive frequencies is returned 'due to symmetry'. Press, et al (2002) make mention of this situation. Be warned that one occasionally sees PSDs defined without this factor two. These, strictly speaking, are called two-sided power spectral densities, but some books are not careful about stating whether one- or two-sided is to be assumed (2002, p. 504). As a result, I infer that spec.pgram() is returning what Press, et al. would call a two-sided PSD. Thus, the power spectrum returned by spec.pgram() should be multiplied by 2 between (DC, nyquist) non-inclusive, and the mean can be resolved separately as the DC component is not returned by spec.pgram(), as: ## N/2 used here because the length of PSD$spec is N/2 rather than N/2 + 1 # as it does not return a DC value. mean(a) + PSD$spec[floor(N/2)] + 2 * sum(PSD$spec[2:floor(N/2)-1]) [1] 64.54347 This situation is closer, and indicates that the normalization is likely too the sum squared amplitude of a[], but it is not EXACT, as I would expect. But I also checked a manual PSD estimate just to show the power spectrum of a single periodicity would actually match the mean or sum squared amplitude of a[] (the discrete function) exactly. a.fft - fft(a) # Two-sided estimate normed to SS Amplitude 1/N * sum(Mod(a.fft[1:floor(N/2)+1])^2) [1] 32 # One-sided estimate normed to SS Amplitude. The first part gets the # quantities across the nyquist range from 0 too fc. The next 2 lines # take out the factor 2 from DC and nyquist since those 2 terms are not # doubled in the one-sided estimate. a.PSD - 1/N * 2 * Mod(a.fft[1:(floor(N/2)+1)])^2 a.PSD[1] - a.PSD[1]/2 #DC freq a.PSD[floor(N/2)+1]- a.PSD[floor(N/2)+1]/2 sum(a.PSD) [1] 64 This proves that the integral of the one-sided PSD estimate of the discrete function a[] is actually exactly the same as the sum squared amplitude of a[], at least in this simple case. Likewise when the periodogram is normalized to the mean-squared amplitude of a[], sum(a.PSD)/N # MS amplitude of a[] [1] 0.5 ## So I don't have to dimension a.PSD[], I took twice the modulus squared # of f[0] too f[c] inclusive all at once, and took out the erroneous # factor 2 from f[1] and f[c] independently a.PSD - 1/(N^2) * 2 * Mod(a.fft[1:(floor(N/2)+1)])^2 a.PSD[1] - a.PSD[1]/2 a.PSD[floor(N/2)+1]- a.PSD[floor(N/2)+1]/2 Of note, is that I have yet to encounter a case where the normalized raw periodogram didn't equal the quantity it was normalized too, even when segmenting data into multiple PSD estimates, and taking their average at each frequency. I have not applied a 'window function' (Welch, Hanning, etc) yet, but the only case I encountered where it's not exact is with overlapping segments, and that, I suspect, is because the first section of the first window, and last section of the last window are only used once as opposed to twice like all of the other data points. If the last segment wrapped around too the first in order too form the last window, I'd expect the norming to be exact again. I have yet to encounter an estimate returned by
Re: [R] R-help Digest, Vol 35, Issue 24
Dear Prof Ripley, First of all, unless you are an english professor, then I do not think you have any business policing language. I'm still very much a student, both in R, and regarding signal analysis. My competence on the subject as compared too your own level of expertise, or my spelling for that matter, may be a contension for you, but it would have been better had you kept that opinion too yourself. There are plenty of other reasons besides laziness or carelessness that people will consistently error in language use, such as learning disorders, head injuries, and/or vertigo. On the contrary, I am aware of the definition of a periodogram, and I know what the unnormalized periodogram in the data I presented looks like. Spec.pgram() is actually normalized too something, because it's discrete integral is not well above the SS amplitude of the signal it computed the periodogram for. In other words, the powers are not in units of around 4,000, which the peak would be if the units were merely the modulus squared of the Fourier coeficients of the data I presented. Alas, the modulus squared of the Fourier coeficients IS the TWO SIDED unnormalized periodogram, ranging from [-fc, fc] | fc=nyquist critical frequency. The definition of the ONE SIDED periodogram IS the modulus squared of the Fourier coeficients ranging over [0, fc], but since the function is even, data points in (0, fc) non-inclusive, need to be multiplied by 2. Thus is according too the definition given by Press, et al (1988, 1992, 2002, c.f. cp 12 13). I'm assuming that R returns an FFT in the same layout as Press, et al describe. Press, et al. are also very clear about the existence of far too many ways of normalizing the periodogram too document, which they stated before delving into particularly how they normalized to the mean squared amplitude of the signal that the periodogram was computed from. In the page before, and perhaps this is where some of the confusion arises from, they document the calculations for MS and SS amplitudes and time integral squared amplitude of the signal in the time domain, not the frequency domain. The page after that, their example only shows how to normalize a periodogram so its sum is equal too the MS amplitude. In short, but starting from SS amplitude: a). sum(a[index=(1:N) or t=(0:N-1)]^2) = SS amplitude calculated in time domain b). 1/N * sum(Mod(fft[-fc:fc])^2) = two sided periodogram that sums too the SS amplitude c). Same as b but over the range [0, fc], and (0, fc) multiplied by 2 is the one sided periodogram, also sums too the SS amplitude For MS amplitude, the procedures are identical, only the time domain is divided by N, and the frequency domain figures are divided by N^2 instead of N. When the periodogram is in power per unit time, as in the above, so that the power is interpretable at N/2+1 independent frequencies, it is a normalized periodogram. spec.pgram() IS normalized, I just do not know what it's normalized too because I can not seem to get spec.pgram to stop tapering (at which point the normalization should be dead on, not just close). By the way, normalized does not automatically mean anything unless to what is stated. I could normalize something arbitrarily to the number of tics on my dogs back side, and still call it normed, or erroneously refer too it as unnormed. If normalized is suposed to mean something specific, then I am confident that more than 90% of undergraduates are not familiar with what the term should mean. Stats and coding and using programs are a human endeavor. This human seems to have made meaning out of terms differently than what those who wrote the documentation seem to have intended. Only, I do not know where the documentation or my understanding may have been missled (R docs, Numerical Recipes, or any other source I looked at since I started). Cheers, KeithC. First, please look up `too' in your dictionary. Second, please study the references on the help page, which give the details. That is what references are for! The references will also answer your question about the reference distribution. The help page does not say it is `normalized' at all: it says it computes the peridogram, and you seem unaware of the definitions of the latter (and beware, there are more than one). On Tue, 24 Jan 2006, Keith Chamberlain wrote: __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] Obtaining the adjusted r-square given the regression
Hi people, I want to obtain the adjusted r-square given a set of coefficients (without the intercept), and I don't know if there is a function that does it. Exist Dear Alexandra, Without knowing what routine you were using that returned an R-Square value too you, it is a little difficult to tell. I am not able to comment on your routine, but I'll try to give a satisfactory response to your first question. Assuming that you used lm() for your linear regression, then all you need to do is see the help page for the lm() summary method, which explains how to get at the values in the summary. In your console type the following: ?summary.lm That help page describes how to get at the values in the summary, and provides an example with getting the full correlation value (whereas the one printed is in a 'pretty' form). I assume the other ones can be accessed similarly, including the adjusted R^2. Your code might look something like the following: # Where object is the result of a call to lm() summary(object)$adj.r.squared The help also has a demonstration for how to exclude the intercept, where you simply add a -1 term too the end of your formula. There may still be an accessor function available, which I do not know about (BUT am certainly interested in if anyone has feedback). In the future, it would help people on the list decode your function when it too is printed in a readable form. The layout of your function in your message may have become garbled when it was stripped of all of its HTML information. If the code was formatted in a readable way when you sent the message, then that's what happened. If that is the case try composing your message in plain text to begin with. I hope this helps. Rgds, KeithC. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] extend.series not zero padding
Dear List, I was trying to verify that I could use extend.series in the wavelets package and kept getting an error when trying to use method=zero. I'm not seeing where my syntax has gone awry. According to the documentation, [see ?extend.series] method: A character string indicating which extension method to use. Possible values are 'periodic', 'reflection', 'zero', 'mean', and 'reflection.inverse'. c-cbind(0:60, 60:0) # setup a series, length will be 61 length(c) [1] 122 dew-extend.series(c,method=zero,length=powerof2, j=log(length(c),2)%/%1) Error in extend.series(c, method = zero, length = powerof2, j = log(length(c), : Invalid argument value for 'method' Other methods work great, such as method=mean. dew-extend.series(c,method=mean,length=powerof2, j=log(length(c),2)%/%1) length(dew) [1] 128 Has this come up in anyone else's experience? If so, what's the workaround so that I can use zero as a method? Rgds, KeithC. __ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
[R] LME data with complicated random correlational structures
Dear List, This is my first post, and I'm a relatively new R user trying to work out a mixed effects model using lme() with random effects, and a correlation structure, and have looked over the archives, R help on lme, corClasses, etc extensively for clues. My programming experience is minimal (1 semester of C). My mentor, who has much more programming experience, but a comparable level of knowledge of R, has been stumped. What follows are 3 questions pertaining to an lme() model, one on the nested hierarcy, 1 on a strategy for a piecewise approach to the variance given I have ~24 hours of data (sampled at 32Hz, 1hr per subject), and one on the corStruct or how to get rid of serial dependencies before lme(). I'm analyzing skin temperature continuously recorded at 32Hz in Baseline (10 min), Testing (~5 min), and Recovery (20 min) epochs of a face recognition experiment. Stimuli are the same in Baseline and Recovery (portrait or landscape), and in testing, participants were tested on their recognition of a list bw portraits presented just before testing started. On some of the portraits 'learned' the eyes were masked, and in others, the eyes were visible. In testing, the portraits have no masking but the stimuli in testing are labeled Eyes and NoEyes. The data structure looks as follows: Subj/Epoch/Stimuli/Time/Temperature There are 8 subjects 9 epochs - 6 of which were just instruction blocks, and one Learning period. Wrt lme(), I figured out how to use subset too isolate just the Baseline, Learning, and Testing Epochs (and avoid epochs with only 1 stimulus level, such as instruction). Data within each epoch are balanced wrt # trials, but not between epochs. Recovery has twice as many trials as Baseline, and Testing has about half. Time for each epoch is roughly that ratio too, although time in each trial differs. Stimuli are the same in Baseline Recovery, but different in Testing, although there are 2 levels in each used epoch. Time Temperature make up the time series, and Temperature is the dependent variable too stimulus. 1- are fixed effects and random effects discrete? That is, if I set something in my model formula as a fixed effect, then it does not make sense to set it as a random effect as well? The documentation (and posts) were not really clear on that point (not that the documentation technically 'should' be per say, just that I got confused). The nested hierarchy for what actually gets analyzed looks as follows: Subj/Epoch/Stimulus/Temperature Reasoning: there are several temperature samples recorded in each trial of Stimulus. Several stimuli in each Epoch, and all of the Epochs for one subject. Subject is random (theoretically) because of sampling in a population, Epoch would be fixed because all participants went through the same sequence of Epochs, but Stimulus varied randomly within an Epoch, which seems inconsistent when I apply it to the lme model as both a fixed and random effect. Temperature ~ Stimulus-1, random=Subj|Subj/Epoch/Stimulus Subset= Epoch==Baseline | Epoch==Testing | Epoch==Recovery I'm looking to correctly allocate error terms for between subjects (Subj) variability, and further delineate the within subject error between Epoch and Stimulus. The current model that I got to work (memory issues namely) is Temperature ~ Stimulus-1, random=Subj|Subj, which I decided to use to get the residuals to have the Subject variability accounted for and subtracted. Would a list of random structures work better? If so, is each item in the list structured just as the random formula? I haven't actually seen/found any examples of a list of random/nesting structures. 2- Is it possible to take a piecewise approach wrt the variance using lme(), such as modeling the variability of each subject first, then using further-nested terms in a model and the residuals from the previous? If so, what caveats exist for interpreting the variances? I'm not interpreting p-values at this point because of another issue. When I try to set up the correlation structure, I run out of memory fast. I've tried this on a mac G5, an HP Pavilion dv1000 (= Pentium 2.6GHz), and a Gateway with an AMD athalon 900MHz processors. Each system has 386M memory or more, one of which has 1G. 3- Is there a way to get rid of the serial dependency BEFORE running the model with LME(), such as initiating a corStruct before placing it in the model? I'm working with so much data that I'm fine with doing the process piecewise. An AR process was difficult because the residuals are not the same length as the data file that I started with. Serial dependencies still gota go, whether via the correlation term in lme() or some other method, because I'll soon be breaking up the variance into components via spectrum(). So I might as well add a 4th. What's the term that gets me too data after AR() has done it's work? I'm thinking that resid() wasn't right but data that the data differ from their original length prior to an AR
