tag:blogger.com,1999:blog-19016597707491708502024-03-05T02:05:27.895-05:00Quantitative EcologyFunctions and tips for use in the R statistical environment.Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.comBlogger34125tag:blogger.com,1999:blog-1901659770749170850.post-80691225195621667202016-12-13T16:04:00.000-05:002016-12-13T16:04:34.472-05:00Recursively search for text in R scripts.Often I find myself needing to search for .R or or .Rmd files in which I have used a specific function. There are a variety of searches that you can do; however, I wanted something that would work recursively at the command line of a Linux or Mac terminal. I found the answer <a href="http://stackoverflow.com/a/28186178">here</a> and have modified it a bit so that it works as a function in Bash (in this example I search through .R and . Rmd files).<br />
<br />
Add the following to your ~/.bash_aliases file and then source ~/.bashrc:<br />
<br />
<div style="overflow: auto;">
<pre></pre>
<pre>findinRfun() {
if [ -z "$1" ]
then
echo "No argument supplied"
else
#search text is passed as argument $1
find . -type f \( -name "*.R" -o -name "*.Rmd" -o -name "*.r" -o -name "*.rmd" \) -print0 | xargs --null grep --with-filename --line-number --no-messages --color --ignore-case "$1"
fi
}
alias findinR=findinRfun
</pre>
<br />
Now, at the command line, you can search through all the files below the current location that contain your search text.<br />
<br />
<pre>findinR mysearchtext
</pre>
<span style="font-size: 85%;"><span style="font-family: "courier new";"></span></span></div>
Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com0tag:blogger.com,1999:blog-1901659770749170850.post-89874018469602392932014-05-16T15:27:00.001-04:002014-05-16T16:22:51.610-04:00Sample uniformly within a fixed radius.I was asked how to do this today and thought that I would share the answer:
<br />
<div style="overflow: auto;">
<pre style='color:#1f1c1b;background-color:#ffffff;'>
<b><span style='color:#b00000;'>## Sample points uniformly within a fixed radius</span></b>
nrand<b><span style='color:#336366;'>=</span></b><span style='color:#b08000;'>1000</span>
maxstep<b><span style='color:#336366;'>=</span></b><span style='color:#b08000;'>10</span>
<b><span style='color:#b00000;'>## Sample data </span></b>
<b><span style='color:#b00000;'>## NB: To get a truly uniform sample over the circle, you must </span></b>
<b><span style='color:#b00000;'>## sample the square of the distance and then transform back.</span></b>
tempdat<b><span style='color:#336366;'><-</span></b><b>data.frame</b>(<span style='color:#0057ae;'>X0=</span><span style='color:#b08000;'>0</span>,<span style='color:#0057ae;'>Y0=</span><span style='color:#b08000;'>0</span>, <span style='color:#0057ae;'>bearing0=</span><span style='color:#b08000;'>0</span>,
<span style='color:#0057ae;'>bad.dist=</span> <b>runif</b>(nrand)<span style='color:#803f00;'>*</span>maxstep,
<span style='color:#0057ae;'>dist2=</span><b>sqrt</b>(<b>runif</b>(nrand)<span style='color:#803f00;'>*</span>maxstep<span style='color:#803f00;'>^</span><span style='color:#b08000;'>2</span>),
<span style='color:#0057ae;'>turningangle=</span><b>runif</b>(nrand)<span style='color:#803f00;'>*</span><span style='color:#b08000;'>2</span><span style='color:#803f00;'>*</span>pi<span style='color:#803f00;'>-</span>pi)
<b><span style='color:#b00000;'>##convert Turning angle to bearing (in this case no change)</span></b>
tempdat<span style='color:#803f00;'>$</span>bearing<b><span style='color:#336366;'>=</span></b>tempdat<span style='color:#803f00;'>$</span>bearing0<span style='color:#803f00;'>+</span>tempdat<span style='color:#803f00;'>$</span>turningangle
<b><span style='color:#b00000;'>## Convert from polar to cartesian coordinates</span></b>
tempdat<span style='color:#803f00;'>$</span>X<b><span style='color:#336366;'><-</span></b>tempdat<span style='color:#803f00;'>$</span>X0<span style='color:#803f00;'>+</span>tempdat<span style='color:#803f00;'>$</span>dist2<span style='color:#803f00;'>*</span><b>sin</b>(tempdat<span style='color:#803f00;'>$</span>bearing)
tempdat<span style='color:#803f00;'>$</span>Y<b><span style='color:#336366;'><-</span></b>tempdat<span style='color:#803f00;'>$</span>Y0<span style='color:#803f00;'>+</span>tempdat<span style='color:#803f00;'>$</span>dist2<span style='color:#803f00;'>*</span><b>cos</b>(tempdat<span style='color:#803f00;'>$</span>bearing)
tempdat<span style='color:#803f00;'>$</span>Xbad<b><span style='color:#336366;'><-</span></b>tempdat<span style='color:#803f00;'>$</span>X0<span style='color:#803f00;'>+</span>tempdat<span style='color:#803f00;'>$</span>bad.dist<span style='color:#803f00;'>*</span><b>sin</b>(tempdat<span style='color:#803f00;'>$</span>bearing)
tempdat<span style='color:#803f00;'>$</span>Ybad<b><span style='color:#336366;'><-</span></b>tempdat<span style='color:#803f00;'>$</span>Y0<span style='color:#803f00;'>+</span>tempdat<span style='color:#803f00;'>$</span>bad.dist<span style='color:#803f00;'>*</span><b>cos</b>(tempdat<span style='color:#803f00;'>$</span>bearing)
<b><span style='color:#b00000;'>##make plots</span></b>
<b>png</b>(<span style='color:#0057ae;'>filename=</span><span style='color:#bf0303;'>"sampleplots.png"</span>,<span style='color:#0057ae;'>width=</span><span style='color:#b08000;'>500</span>,<span style='color:#0057ae;'>height=</span><span style='color:#b08000;'>1000</span>)
<b>par</b>(<span style='color:#0057ae;'>mfrow=</span><b>c</b>(<span style='color:#b08000;'>2</span>,<span style='color:#b08000;'>1</span>))
<b>plot</b>(Ybad<span style='color:#803f00;'>~</span>Xbad, <span style='color:#0057ae;'>data=</span>tempdat, <span style='color:#0057ae;'>asp=</span><span style='color:#b08000;'>1</span>, <span style='color:#0057ae;'>main=</span><span style='color:#bf0303;'>"Center is oversampled"</span>)
<b>plot</b>(Y<span style='color:#803f00;'>~</span>X, <span style='color:#0057ae;'>data=</span>tempdat, <span style='color:#0057ae;'>asp=</span><span style='color:#b08000;'>1</span>, <span style='color:#0057ae;'>main=</span><span style='color:#bf0303;'>"Uniform across space"</span>)
<b>dev.off</b>()
</pre>
</pre>
<span style="font-size: 85%;"><span style="font-family: courier new;"></span></span></div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiU7ZvNH49RoPrQkPB6jkzVey-QQSOhgiNcqhtwZEHMq8Vzru2BH7UsE6pu57eS4F1drhXfz8JWmRnQjHwiP3Fi0DIbHL11BNW0k8h_HSVijgO2fKMk5c_aNtk1x62Y_t9bgsgcMLFy5hOP/s1600/sampleplots.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiU7ZvNH49RoPrQkPB6jkzVey-QQSOhgiNcqhtwZEHMq8Vzru2BH7UsE6pu57eS4F1drhXfz8JWmRnQjHwiP3Fi0DIbHL11BNW0k8h_HSVijgO2fKMk5c_aNtk1x62Y_t9bgsgcMLFy5hOP/s1600/sampleplots.png" height="640" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
</div>
Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com1tag:blogger.com,1999:blog-1901659770749170850.post-34477763536005476122014-04-28T17:03:00.000-04:002014-10-20T13:14:19.042-04:00RStudio and X2GoAfter a recent system upgrade (to <a href="http://www.linuxmint.com/download_lmde.php" target="_blank">Linux Mint LMDE</a>), I was no longer able to run <a href="https://www.rstudio.com/" target="_blank">RStudio</a> through the remote desktop application <a href="http://wiki.x2go.org/doku.php" target="_blank">X2Go</a>. It turns out that this is due to a problem with the Qt libraries (see this <a href="https://hpc2.org/kb/development/deploying-qt-applications">website</a>).
As suggested <a href="https://forum.kde.org/viewtopic.php?f=66&t=110934">here</a>, I just deleted all the Qt libraries used by RStudio:<br />
<br />
<div style="overflow: auto;">
<pre>sudo rm /usr/lib/rstudio/bin/libQt* </pre>
<pre> </pre>
<span style="font-size: 85%;"><span style="font-family: courier new;"></span></span></div>
<br />
and all was well. (NB, the specific link to these files will vary depending on the flavor of Linux you use.)<br />
<br />
**EDIT: If you are installing this on pure Debian (I just confirmed this with a new Jessie install) you will need to install appropriate libraries too:<br />
<br />
<div style="overflow: auto;">
<pre>sudo aptitude install libqtwebkit4</pre>
<pre> </pre>
</div>
<br />
<br />
<br />
<br />
Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com0tag:blogger.com,1999:blog-1901659770749170850.post-16292900287251715352010-11-10T15:56:00.002-05:002010-11-10T16:00:39.429-05:00At last... I have been suffering with XEmacs displaying odd characters instead of the quotation marks that are used in R help files. This was driving me up the wall because it makes the files (and R output in general) very hard to read; however, I finally diagnosed the problem: Xemacs was not recognizing UTF-8 encoding. Below is a quote from <a href="http://www.maruko.ca/i18n/">Marjan Parsa</a> that describes how to set up Emacs and XEmacs to automatically detect UTF-8 files. My quality of life has already improved.<br /><br /><br /><blockquote>How can I get XEmacs to work with UTF-8 files?<br /><br /> * Set up XEmacs so that it autodetects UTF-8 encoded files.<br /> * In the case of starting a new file in a non-UTF-8 locale, set the file coding system to UTF-8 using C-x RET f.<br /> * If running XEmacs in non-graphical mode in a UTF-8 xterm, set the terminal coding system to UTF-8 using C-x RET t.<br /><br />If you want XEmacs to load UTF-8 files correctly, add the following lines to your ~/.xemacs/init.el:<br /><br />(require 'un-define)<br />(set-coding-priority-list '(utf-8))<br />(set-coding-category-system 'utf-8 'utf-8)<br /><br />Note that Emacs does not deal well with these additions, so if you also run Emacs, then adding the following will keep Emacs from complaining:<br /><br />;; Are we running XEmacs or Emacs?<br />(defvar running-xemacs (string-match "XEmacs\\|Lucid" emacs-version))<br /><br />...<br /><br />(if (not running-xemacs) nil<br /> ;; enable Mule-UCS<br /> (require 'un-define)<br /><br /> ;; by default xemacs does not autodetect Unicode<br /> (set-coding-priority-list '(utf-8))<br /> (set-coding-category-system 'utf-8 'utf-8))<br /><br />These lines will get XEmacs to load UTF-8 files in UTF-8 mode (it will display a "u" in the bottom left corner of your status bar). If you have already loaded a file and would like to start inputting UTF-8, you can use C-x RET f, to set the file coding system to UTF-8. Note that you may additionally have to set the terminal coding system to UTF-8. This seems to be necessary, for example, in the case where XEmacs is run in non-graphical mode inside a UTF-8 enabled xterm. You can set the terminal encoding using C-x RET t.<br /><br />Caution: I have had problems with XEmacs double encoding in the case where 1) the file contains UTF-8, 2) the file is loaded in non-UTF-8 mode, 3) the user switches to UTF-8 mode (using C-x RET f), 4) enters some text, and 5) saves. In other words, if your file already contains UTF-8 characters, make sure that it is loaded in UTF-8 mode before editing it.</blockquote>Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com0tag:blogger.com,1999:blog-1901659770749170850.post-19808833214323947132009-10-22T14:38:00.008-04:002010-02-19T11:29:09.847-05:00ISO weekI am working with a model that produces estimates of snow water equivalent through time. Because I deal with large spatial extents, I decided to have the model produce weekly averages. The problem with this is knowing which file to access for a given date. The snow files are saved using an <a href="http://en.wikipedia.org/wiki/ISO_week_date">ISO-8061 week number</a> (I had no idea how many date "standards" there were, but this is a common one). I thought that getting the week number out of a date format in R would be easy, but I was wrong. It turns out to be platform specific and regardless, as far as I could tell, there is no built-in way to calculate the ISO week. Thus the following function.<br /><br /><div style="overflow: auto;"><pre><span style="color: rgb(20, 19, 18);">ISOweek</span><span style="color: rgb(51, 99, 102);"><b><-</b></span><span style="color: rgb(0, 0, 191);">function</span><span style="color: rgb(20, 19, 18);">(date,</span><span style="color: rgb(0, 87, 174);">format=</span><span style="color: rgb(191, 3, 3);">"%Y-%m-%d"</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(0, 87, 174);">return.val=</span><span style="color: rgb(191, 3, 3);">"weekofyear"</span><span style="color: rgb(20, 19, 18);">){</span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##converts dates into "dayofyear" or "weekofyear", the latter providing the ISO-8601 week</b></span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##date should be a vector of class Date or a vector of formatted character strings</b></span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##format refers to the date form used if a vector of</b></span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>## character strings is supplied</b></span><br /><span style="color: rgb(20, 19, 18);"> </span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##convert date to POSIXt format </b></span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(0, 0, 191);">if</span><span style="color: rgb(20, 19, 18);">(<b>class</b>(date)[</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(20, 19, 18);">]</span><span style="color: rgb(128, 63, 0);">%in%</span><span style="color: rgb(20, 19, 18);"><b>c</b>(</span><span style="color: rgb(191, 3, 3);">"Date"</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(191, 3, 3);">"character"</span><span style="color: rgb(20, 19, 18);">)){</span><br /><span style="color: rgb(20, 19, 18);"> date</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);"><b>as.POSIXlt</b>(date,</span><span style="color: rgb(0, 87, 174);">format=</span><span style="color: rgb(20, 19, 18);">format)</span><br /><span style="color: rgb(20, 19, 18);"> }</span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(0, 0, 191);">if</span><span style="color: rgb(20, 19, 18);">(<b>class</b>(date)[</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(20, 19, 18);">]</span><span style="color: rgb(128, 63, 0);">!=</span><span style="color: rgb(191, 3, 3);">"POSIXt"</span><span style="color: rgb(20, 19, 18);">){</span><br /><span style="color: rgb(20, 19, 18);"> <b>print</b>(</span><span style="color: rgb(191, 3, 3);">"Date is of wrong format."</span><span style="color: rgb(20, 19, 18);">)</span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(0, 0, 191);">break</span><br /><span style="color: rgb(20, 19, 18);"> }</span>else if(class(date)[2]=="POSIXct"){<br /> date=as.POSIXlt(date)<br />}<br /><br /><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(0, 0, 191);">if</span><span style="color: rgb(20, 19, 18);">(return.val</span><span style="color: rgb(128, 63, 0);">==</span><span style="color: rgb(191, 3, 3);">"dayofyear"</span><span style="color: rgb(20, 19, 18);">){</span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##add 1 because POSIXt is base zero</b></span><br /><span style="color: rgb(20, 19, 18);"> <b>return</b>(date</span><span style="color: rgb(128, 63, 0);">$</span><span style="color: rgb(20, 19, 18);">yday</span><span style="color: rgb(128, 63, 0);">+</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(20, 19, 18);">)</span><br /><span style="color: rgb(20, 19, 18);"> }</span><span style="color: rgb(0, 0, 191);">else</span><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(0, 0, 191);">if</span><span style="color: rgb(20, 19, 18);">(return.val</span><span style="color: rgb(128, 63, 0);">==</span><span style="color: rgb(191, 3, 3);">"weekofyear"</span><span style="color: rgb(20, 19, 18);">){</span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##Based on the ISO8601 weekdate system,</b></span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>## Monday is the first day of the week</b></span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>## W01 is the week with 4 Jan in it.</b></span><br /><span style="color: rgb(20, 19, 18);"> year</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(176, 128, 0);">1900</span><span style="color: rgb(128, 63, 0);">+</span><span style="color: rgb(20, 19, 18);">date</span><span style="color: rgb(128, 63, 0);">$</span><span style="color: rgb(20, 19, 18);">year</span><br /><span style="color: rgb(20, 19, 18);"> jan4</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);"><b>strptime</b>(<b>paste</b>(year,</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(176, 128, 0);">4</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(0, 87, 174);">sep=</span><span style="color: rgb(191, 3, 3);">"-"</span><span style="color: rgb(20, 19, 18);">),</span><span style="color: rgb(0, 87, 174);">format=</span><span style="color: rgb(191, 3, 3);">"%Y-%m-%d"</span><span style="color: rgb(20, 19, 18);">)</span><br /><span style="color: rgb(20, 19, 18);"> wday</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">jan4</span><span style="color: rgb(128, 63, 0);">$</span><span style="color: rgb(20, 19, 18);">wday</span><br /><span style="color: rgb(20, 19, 18);"> </span><br /><span style="color: rgb(20, 19, 18);"> wday[wday</span><span style="color: rgb(128, 63, 0);">==</span><span style="color: rgb(176, 128, 0);">0</span><span style="color: rgb(20, 19, 18);">]</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(176, 128, 0);">7</span><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##convert to base 1, where Monday == 1, Sunday==7</b></span><br /><span style="color: rgb(20, 19, 18);"> </span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##calculate the date of the first week of the year</b></span><br /><span style="color: rgb(20, 19, 18);"> weekstart</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">jan4</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);">(wday</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(20, 19, 18);">)</span><span style="color: rgb(128, 63, 0);">*</span><span style="color: rgb(176, 128, 0);">86400</span><span style="color: rgb(20, 19, 18);"> </span><br /><span style="color: rgb(20, 19, 18);"> weeknum</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);"><b>ceiling</b>(<b>as.numeric</b>((<b>difftime</b>(date,weekstart,</span><span style="color: rgb(0, 87, 174);">units=</span><span style="color: rgb(191, 3, 3);">"days"</span><span style="color: rgb(20, 19, 18);">)</span><span style="color: rgb(128, 63, 0);">+</span><span style="color: rgb(176, 128, 0);">0.1</span><span style="color: rgb(20, 19, 18);">)</span><span style="color: rgb(128, 63, 0);">/</span><span style="color: rgb(176, 128, 0);">7</span><span style="color: rgb(20, 19, 18);">))</span><br /><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>#########################################################################</b></span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##calculate week for days of the year occuring in the next year's week 1.</b></span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>#########################################################################</b></span><br /><span style="color: rgb(20, 19, 18);"> mday</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">date</span><span style="color: rgb(128, 63, 0);">$</span><span style="color: rgb(20, 19, 18);">mday</span><br /><span style="color: rgb(20, 19, 18);"> wday</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">date</span><span style="color: rgb(128, 63, 0);">$</span><span style="color: rgb(20, 19, 18);">wday</span><br /><span style="color: rgb(20, 19, 18);"> wday[wday</span><span style="color: rgb(128, 63, 0);">==</span><span style="color: rgb(176, 128, 0);">0</span><span style="color: rgb(20, 19, 18);">]</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(176, 128, 0);">7</span><br /><span style="color: rgb(20, 19, 18);"> year</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);"><b>ifelse</b>(weeknum</span><span style="color: rgb(128, 63, 0);">==</span><span style="color: rgb(176, 128, 0);">53</span><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(128, 63, 0);">&</span><span style="color: rgb(20, 19, 18);"> mday</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);">wday</span><span style="color: rgb(128, 63, 0);">>=</span><span style="color: rgb(176, 128, 0);">28</span><span style="color: rgb(20, 19, 18);">,year</span><span style="color: rgb(128, 63, 0);">+</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(20, 19, 18);">,year)</span><br /><span style="color: rgb(20, 19, 18);"> weeknum</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);"><b>ifelse</b>(weeknum</span><span style="color: rgb(128, 63, 0);">==</span><span style="color: rgb(176, 128, 0);">53</span><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(128, 63, 0);">&</span><span style="color: rgb(20, 19, 18);"> mday</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);">wday</span><span style="color: rgb(128, 63, 0);">>=</span><span style="color: rgb(176, 128, 0);">28</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(20, 19, 18);">,weeknum)</span><br /><span style="color: rgb(20, 19, 18);"> </span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>################################################################</b></span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##calculate week for days of the year occuring prior to week 1.</b></span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>################################################################</b></span><br /><span style="color: rgb(20, 19, 18);"> </span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##first calculate the numbe of weeks in the previous year</b></span><br /><span style="color: rgb(20, 19, 18);"> year.shift</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">year</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(176, 128, 0);">1</span><br /><span style="color: rgb(20, 19, 18);"> jan4.shift</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);"><b>strptime</b>(<b>paste</b>(year.shift,</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(176, 128, 0);">4</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(0, 87, 174);">sep=</span><span style="color: rgb(191, 3, 3);">"-"</span><span style="color: rgb(20, 19, 18);">),</span><span style="color: rgb(0, 87, 174);">format=</span><span style="color: rgb(191, 3, 3);">"%Y-%m-%d"</span><span style="color: rgb(20, 19, 18);">)</span><br /><span style="color: rgb(20, 19, 18);"> wday</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">jan4.shift</span><span style="color: rgb(128, 63, 0);">$</span><span style="color: rgb(20, 19, 18);">wday</span><br /><span style="color: rgb(20, 19, 18);"> wday[wday</span><span style="color: rgb(128, 63, 0);">==</span><span style="color: rgb(176, 128, 0);">0</span><span style="color: rgb(20, 19, 18);">]</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(176, 128, 0);">7</span><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##convert to base 1, where Monday == 1, Sunday==7</b></span><br /><span style="color: rgb(20, 19, 18);"> weekstart</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">jan4.shift</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);">(wday</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(20, 19, 18);">)</span><span style="color: rgb(128, 63, 0);">*</span><span style="color: rgb(176, 128, 0);">86400</span><br /><span style="color: rgb(20, 19, 18);"> weeknum.shift</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);"><b>ceiling</b>(<b>as.numeric</b>((<b>difftime</b>(date,weekstart)</span><span style="color: rgb(128, 63, 0);">+</span><span style="color: rgb(176, 128, 0);">0.1</span><span style="color: rgb(20, 19, 18);">)</span><span style="color: rgb(128, 63, 0);">/</span><span style="color: rgb(176, 128, 0);">7</span><span style="color: rgb(20, 19, 18);">))</span><br /><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##update year and week</b></span><br /><span style="color: rgb(20, 19, 18);"> year</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);"><b>ifelse</b>(weeknum</span><span style="color: rgb(128, 63, 0);">==</span><span style="color: rgb(176, 128, 0);">0</span><span style="color: rgb(20, 19, 18);">,year.shift,year)</span><br /><span style="color: rgb(20, 19, 18);"> weeknum</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);"><b>ifelse</b>(weeknum</span><span style="color: rgb(128, 63, 0);">==</span><span style="color: rgb(176, 128, 0);">0</span><span style="color: rgb(20, 19, 18);">,weeknum.shift,weeknum)</span><br /><span style="color: rgb(20, 19, 18);"> </span><br /><span style="color: rgb(20, 19, 18);"> <b>return</b>(<b>list</b>(</span><span style="color: rgb(191, 3, 3);">"year"</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">year,</span><span style="color: rgb(191, 3, 3);">"weeknum"</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">weeknum))</span><br /><span style="color: rgb(20, 19, 18);"> }</span><span style="color: rgb(0, 0, 191);">else</span><span style="color: rgb(20, 19, 18);">{</span><br /><span style="color: rgb(20, 19, 18);"> <b>print</b>(</span><span style="color: rgb(191, 3, 3);">"Unknown return.val"</span><span style="color: rgb(20, 19, 18);">)</span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(0, 0, 191);">break</span><br /><span style="color: rgb(20, 19, 18);"> }</span><br /><span style="color: rgb(20, 19, 18);">}</span><br /><br /><br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br /><br />Of course after I wrote this function, I found <a href="http://tinyurl.com/yj2l6rc">this</a> thread in R-help. Gustaf Rydevik provided a function that is similar to mine (his function is also about twice as fast); however, on my computer it was giving incorrect years and week numbers for days in January that occur in the previous year's final week (e.g., 1 January 2010 should be week 53 of 2009). Reading through the rest of the thread indicates that the confusion involving week numbers is widespread. I guess the moral of the story is pick a standard and stick to it.<br /><br />Please let me know in the comments if you encounter any problems running this function (or if you can suggest ways to make it more efficient).Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com7tag:blogger.com,1999:blog-1901659770749170850.post-84641183471860349202009-10-22T14:18:00.002-04:002009-10-22T14:23:43.056-04:00Leap yearsA quick function that when provided a numeric vector of years returns a boolean vector where TRUE == Leap year.<br /><div style="overflow: auto;"><pre><br /><span style='color: #141312'>is.leapyear</span><span style='color: #336366'><b>=</b></span><span style='color: #0000bf'>function</span><span style='color: #141312'>(year){</span><br /><span style='color: #141312'> </span><span style='color: #888786'><i>#http://en.wikipedia.org/wiki/Leap_year</i></span><br /><span style='color: #141312'> <b>return</b>(((year </span><span style='color: #803f00'>%%</span><span style='color: #141312'> </span><span style='color: #b08000'>4</span><span style='color: #141312'> </span><span style='color: #803f00'>==</span><span style='color: #141312'> </span><span style='color: #b08000'>0</span><span style='color: #141312'>) </span><span style='color: #803f00'>&</span><span style='color: #141312'> (year </span><span style='color: #803f00'>%%</span><span style='color: #141312'> </span><span style='color: #b08000'>100</span><span style='color: #141312'> </span><span style='color: #803f00'>!=</span><span style='color: #141312'> </span><span style='color: #b08000'>0</span><span style='color: #141312'>)) </span><span style='color: #803f00'>|</span><span style='color: #141312'> (year </span><span style='color: #803f00'>%%</span><span style='color: #141312'> </span><span style='color: #b08000'>400</span><span style='color: #141312'> </span><span style='color: #803f00'>==</span><span style='color: #141312'> </span><span style='color: #b08000'>0</span><span style='color: #141312'>))</span><br /><span style='color: #141312'>}</span></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div>Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com1tag:blogger.com,1999:blog-1901659770749170850.post-65394931272873544582009-09-03T08:17:00.006-04:002009-09-03T09:11:07.669-04:00Truncated Normal DistributionMany distributions may be used to describe patterns that are non-negative; however, there are not as many choices when an upper bound is also needed (although the <a href="http://en.wikipedia.org/wiki/Beta_distribution">beta distribution</a> <span style="text-decoration: underline;"></span>is very flexible). For various reasons, truncated distributions are sometimes preferred, and the <a href="http://en.wikipedia.org/wiki/Truncated_normal_distribution">truncated normal</a> is particularly popular. While R has a package that includes the standard functions for this distribution (see rtnorm, dtnorm, etc. in the <a href="http://cran.r-project.org/web/packages/msm/index.html">msm pacakge</a>), the true expectation and variance of the distribution may be of interest. It turns out that the first two moments of the truncated normal are not too hard to calculate (but worth writing functions for):<br /><div style="overflow: auto;"><pre><pre><br /><span style="color: rgb(20, 19, 18);">mean.tnorm</span><span style="color: rgb(51, 99, 102);"><b><-</b></span><span style="color: rgb(0, 0, 191);">function</span><span style="color: rgb(20, 19, 18);">(mu,sd,lower,upper){</span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##return the expectation of a truncated normal distribution</b></span><br /><span style="color: rgb(20, 19, 18);"> lower.std</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">(lower</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);">mu)</span><span style="color: rgb(128, 63, 0);">/</span><span style="color: rgb(20, 19, 18);">sd</span><br /><span style="color: rgb(20, 19, 18);"> upper.std</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">(upper</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);">mu)</span><span style="color: rgb(128, 63, 0);">/</span><span style="color: rgb(20, 19, 18);">sd</span><br /><span style="color: rgb(20, 19, 18);"> mean</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">mu</span><span style="color: rgb(128, 63, 0);">+</span><span style="color: rgb(20, 19, 18);">sd</span><span style="color: rgb(128, 63, 0);">*</span><span style="color: rgb(20, 19, 18);">(<b>dnorm</b>(lower.std)</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);"><b>dnorm</b>(upper.std))</span><span style="color: rgb(128, 63, 0);">/</span><span style="color: rgb(20, 19, 18);"><br /> (<b>pnorm</b>(upper.std)</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);"><b>pnorm</b>(lower.std))</span><br /><span style="color: rgb(20, 19, 18);"> <b>return</b>(mean)</span><br /><span style="color: rgb(20, 19, 18);">}</span><br /><br /><span style="color: rgb(20, 19, 18);">var.tnorm</span><span style="color: rgb(51, 99, 102);"><b><-</b></span><span style="color: rgb(0, 0, 191);">function</span><span style="color: rgb(20, 19, 18);">(mu,sd,lower,upper){</span><br /><span style="color: rgb(20, 19, 18);"> </span><span style="color: rgb(176, 0, 0);"><b>##return the variance of a truncated normal distribution</b></span><br /><span style="color: rgb(20, 19, 18);"> lower.std</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">(lower</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);">mu)</span><span style="color: rgb(128, 63, 0);">/</span><span style="color: rgb(20, 19, 18);">sd</span><br /><span style="color: rgb(20, 19, 18);"> upper.std</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">(upper</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);">mu)</span><span style="color: rgb(128, 63, 0);">/</span><span style="color: rgb(20, 19, 18);">sd</span><br /><span style="color: rgb(20, 19, 18);"> variance</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);">sd</span><span style="color: rgb(128, 63, 0);">^</span><span style="color: rgb(176, 128, 0);">2</span><span style="color: rgb(128, 63, 0);">*</span><span style="color: rgb(20, 19, 18);">(</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(128, 63, 0);">+</span><span style="color: rgb(20, 19, 18);">(lower.std</span><span style="color: rgb(128, 63, 0);">*</span><span style="color: rgb(20, 19, 18);"><b>dnorm</b>(lower.std)</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);">upper.std</span><span style="color: rgb(128, 63, 0);">*</span><span style="color: rgb(20, 19, 18);"><b>dnorm</b>(upper.std))</span><span style="color: rgb(128, 63, 0);">/</span><br /><span style="color: rgb(20, 19, 18);"> (<b>pnorm</b>(upper.std)</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);"><b>pnorm</b>(lower.std))</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);">((<b>dnorm</b>(lower.std)</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);"><b>dnorm</b>(upper.std))</span><span style="color: rgb(128, 63, 0);">/</span><br /><span style="color: rgb(20, 19, 18);"> (<b>pnorm</b>(upper.std)</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(20, 19, 18);"><b>pnorm</b>(lower.std)))</span><span style="color: rgb(128, 63, 0);">^</span><span style="color: rgb(176, 128, 0);">2</span><span style="color: rgb(20, 19, 18);">)</span><br /><span style="color: rgb(20, 19, 18);"> <b>return</b>(variance)</span><br /><span style="color: rgb(20, 19, 18);">}</span><br /><br /><span style="color: rgb(176, 0, 0);"><b>###Testing</b></span><br /><span style="color: rgb(20, 19, 18);">> <b>library</b>(msm)</span><br /><span style="color: rgb(20, 19, 18);">> a</span><span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(20, 19, 18);"><b>rtnorm</b>(</span><span style="color: rgb(176, 128, 0);">1000000</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(176, 128, 0);">5</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(176, 128, 0);">2</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(176, 128, 0);">3</span><span style="color: rgb(20, 19, 18);">)</span><br /><span style="color: rgb(20, 19, 18);">> <b>paste</b>(<b>mean</b>(a),<b>var</b>(a))</span><br /><span style="color: rgb(20, 19, 18);">[</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(20, 19, 18);">] </span><span style="color: rgb(191, 3, 3);">"1.52135857341077 0.197281057170982"</span><br /><span style="color: rgb(20, 19, 18);">> <b>paste</b>(<b>mean.tnorm</b>(</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(176, 128, 0);">5</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(176, 128, 0);">2</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(176, 128, 0);">3</span><span style="color: rgb(20, 19, 18);">),<b>var.tnorm</b>(</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(176, 128, 0);">5</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(176, 128, 0);">2</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(20, 19, 18);">,</span><span style="color: rgb(176, 128, 0);">3</span><span style="color: rgb(20, 19, 18);">))</span><br /><span style="color: rgb(20, 19, 18);">[</span><span style="color: rgb(176, 128, 0);">1</span><span style="color: rgb(20, 19, 18);">] </span><span style="color: rgb(191, 3, 3);">"1.52090857118 0.197111175109889"</span></pre><br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div>Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com1tag:blogger.com,1999:blog-1901659770749170850.post-78861918538469388252009-08-28T16:50:00.006-04:002009-08-28T17:04:01.620-04:00DropboxI have been trying to convince my collaborators to use git so that we all work on the same page; however, I have accepted the fact that not everyone wants to get that involved. For these collaborations, I have started relying on a free web-based service called <a href="https://www.getdropbox.com/referrals/NTExNzc0NDU5">Dropbox</a>. The free version gives you a 2GB account and allows you to securely share files with others without the hassle of setting up accounts on an ftp server. There is even a client for Linux (thank you!).<br /><br />On Linux, after downloading and installing the program, you just have to start the daemon:<div style="overflow: auto;"><pre>$dropbox start -i<br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div>During the set-up process, a folder is created in your home directory called "Dropbox" and everything you put in there is automatically synced with all the other computers you have linked to the account. You can explicitly share folders with others, or place files in the "Public" subdirectory and then share the associated url.<br /><br />I have started using this service as a virtual flash drive and as an alternative to emailing enormous files.Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com0tag:blogger.com,1999:blog-1901659770749170850.post-41483662379908391462009-08-28T12:07:00.006-04:002009-11-26T09:44:50.332-05:00git cheat sheetFrom what I have seen, the best description of what Git is and how it works operationally was written by <a href="http://www.eecs.harvard.edu/%7Ecduan/technical/git/">Charles Duan</a>. Given that Charles has done such a great job, I am just going to summarize a few key commands here (also based on <a href="http://toolmantim.com/articles/setting_up_a_new_remote_git_repository">this</a> tutorial).<br /><br /><span style="font-weight: bold;">A) Create a local git repo out of your current project directory.</span><br /><div style="overflow: auto;"><pre><br />$ cd ~/myproj<br />$ git init<br />$ git add .<br />$ git commit -m "first commit"<br /></pre></div><br /><span style="font-weight: bold;">B) Create an empty remote git repo and push current project repo there.</span> <div style="overflow: auto;"><pre><br />$ ssh my.remote.server<br /><br />> mkdir gitroot<br />> mkdir gitroot/myproj<br />> cd gitroot/myproj<br />> git --bare init<br />> exit<br /><br />$ cd ~/myproj<br />$ git remote add origin ssh://my.remote.server/home/username/gitroot/myproj<br />$ git push origin master<br /></pre></div><br /><span style="font-weight: bold;">C) Check out a remote git repo and send updates back.</span><br /><div style="overflow: auto;"><pre>$ git clone ssh://myserver.com/home/username/gitroot/myproj<br /><make><br />$ git push<br /></make></pre></div><span style="font-weight: bold;">D) Create and switch among branches</span><br /><div style="overflow: auto;"><pre>$ git branch [new-head-name] [reference-to-branching-point]<br /></pre></div>For example:<div style="overflow: auto;"><pre>$ git branch testing HEAD<br /></pre></div>You switch among branches simply by checking them out from your own repo:<br /><div style="overflow: auto;"><pre>$ git checkout testing<br />$ git push origin local-branch-name<br /></pre></div><span style="font-weight: bold;">E) Merge branches</span>.<span style="font-weight: bold;"><span style="font-weight: bold;"></span></span> Again, look <a href="http://www.eecs.harvard.edu/%7Ecduan/technical/git/git-3.shtml">here</a> for a detailed description of how this works in practice.<br /><div style="overflow: auto;"><pre>$ git checkout master<br />$ git merge testing<br /></pre></div><span style="font-weight: bold;">F) Convert svn repos to git. </span>There are <a href="http://pauldowman.com/2008/07/26/how-to-convert-from-subversion-to-git/">other ways to do this</a> if you want to preserve the tree of your svn repo; however, if all you want is to strip out the .svn folders:<br /><div style="overflow: auto;"><pre>$cd ~/mysvnproj<br />$ svn export . ../myproj<br /></pre></div>The problem is that only files under version control are exported. If you just want to remove svn version control from your current directory (i.e., recursively strip out the .svn folders), then use:<br /><pre>$find . -name ".svn" -type d -exec rm -rf {} \;</pre>Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com0tag:blogger.com,1999:blog-1901659770749170850.post-77588307610568363532009-08-19T11:02:00.026-04:002009-08-31T09:59:25.847-04:00SVN gets got by gitMy Subversion usage finally stabilized, but I realized that my major problem was that commits were relatively expensive (i.e., in order for me to make any commit I had to open an ssh connection to my remote server and send in my updates). Generally this is no big deal, but it made me less likely to commit on a regular basis -- especially when at home or travelling -- so I was falling back into the bad habit of making new copies of code with increasingly obtuse names (myscript.R -> myscript2.R -> myscript2-working.R, etc.). I guess I could develop my own personal file nomenclature system, but frankly, that is what version control software is for.<br /><br />Earlier this year I stumbled upon an <a href="http://craphound.com/?p=2171">article</a> by Cory Doctorow in which he describes how he uses a Python script (<a href="http://bitbucketlabs.net/flashbake/">Flashbake</a>) to automate the <a href="http://git-scm.com/">git version control system</a>. This led me to a <a href="http://www.youtube.com/watch?v=4XpnKHJAok8">video</a> of <a href="http://en.wikipedia.org/wiki/Linus_Torvalds">Linus Torvalds</a> preaching the Gospel of Git at the Googleplex. The chief architect of the Linux kernel is known to be fairly, um... opinionated, so I was prepared to take his comments well salted; however, I ultimately found his arguments for git compelling and decided to give git a whirl.<br /><br />Git is many things to many people, but my reasons to boot svn to the curb are that I want to:<br /><br />1) Make local commits often.<br />2) Branch at will.<br />3) Switch among branches easily.<br />4) Merge branches easily.<br />5) Share repos among collaborators.<br /><br />For my workflow, git does all of these things better than svn. While the ability to perform local commits was the killer feature for me, I have come to appreciate branching and merging (something I generally avoided with svn due to the often disastrous outcomes). So there it is. I like Git (and not <span style="font-style: italic;">just</span> because <a href="http://unethicalblogger.com/posts/2009/01/im_using_git_because_it_makes_me_feel_cool">it makes me feel cool!</a>).<br /><br />A friend of mine just mentioned that he likes <a href="http://bazaar-vcs.org/">bazaar version control</a>. It looks very similar to git, (and may have some <a href="http://bazaar-vcs.org/BzrVsGit">advantages</a> for certain situations -- especially if Windows users are loathe to install cygwin); however, I am already committed (for at least a few months).Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com0tag:blogger.com,1999:blog-1901659770749170850.post-62853538654783700152009-07-20T18:40:00.010-04:002009-08-21T10:40:30.636-04:00When subversion does not add up.I have been using Subversion to manage all of my active projects -- it is a fantastic tool but I ran into a bit of a problem recently when I accidentally added several hundred log files as I was preparing for a commit (using the automated method I describe <a href="http://quantitative-ecology.blogspot.com/search/label/subversion">here</a> -- so <span style="font-style: italic;">caveat emptor</span>). Most suggestions I found on the internet were some variation of:<br /><div style="overflow: auto;"><pre>$ svn --recursive revert .<br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div>The problem with this solution is that all uncommitted changes to existing files (i.e., modified files that were already present in the repository) would be lost. Ouch! I discovered <a href="http://svn.haxx.se/users/archive-2008-01/0307.shtml">this</a> solution; however, it did not work for me because I had already made the mistake of deleting the offending files. Finally, I modified the suggested command and all is well:<br /><div style="overflow: auto;"><pre>$ svn st | grep log > badfile.txt<br />$ sed 's/! /.\//' badfile.txt > newfile.txt<br />$ less newfile.txt | xargs svn revert<br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div>Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com0tag:blogger.com,1999:blog-1901659770749170850.post-14346762270654842382008-12-13T16:44:00.008-05:002008-12-13T17:59:54.153-05:00R matrices in C functionsUsing the .C() function in R, you can only pass vectors. Since R stores matrices columnwise as vectors anyhow, they can be passed to your C function as vectors (along with the number of rows in the matrix) and then accessed in familiar [row,col] manner using the following C functions (idea from <a href="https://stat.ethz.ch/pipermail/r-help/2005-April/070488.html">here</a>):<br /><div style="overflow: auto;"><pre><span style="color: rgb(128, 0, 0);">int</span> Cmatrix(<span style="color: rgb(128, 0, 0);">int</span> *row, <span style="color: rgb(128, 0, 0);">int</span> *col, <span style="color: rgb(128, 0, 0);">int</span> *nrows){<br /><span style="color: rgb(128, 128, 128);"><i>/* Converts row-col notation into base-zero vector notation, based on a column-wise conversion*/</i></span><br /><span style="color: rgb(128, 0, 0);">int</span> vector_loc;<br />vector_loc=(*row)-<span style="color: rgb(0, 0, 255);">1</span> + ((*col)-<span style="color: rgb(0, 0, 255);">1</span>)*(*nrows);<br /><b>return</b>(vector_loc);<br />}<br /><br /><span style="color: rgb(128, 0, 0);">void</span> Rmatrix(<span style="color: rgb(128, 0, 0);">int</span> *vector_loc,<span style="color: rgb(128, 0, 0);">int</span> *row, <span style="color: rgb(128, 0, 0);">int</span> *col, <span style="color: rgb(128, 0, 0);">int</span> *nrows){<br /><span style="color: rgb(128, 128, 128);"><i>/*Converts vector notation into row-col notation*/</i></span><br /><span style="color: rgb(128, 128, 128);"><i>/*vector_loc is the vector address of the matrix (base zero)*/</i></span><br /><span style="color: rgb(128, 128, 128);"><i>/*row and col are pointers to the row and col variables that will be updated */</i></span><br />*col=floor(*vector_loc / *nrows)+<span style="color: rgb(0, 0, 255);">1</span>;<br />*row=*vector_loc-((*col)-<span style="color: rgb(0, 0, 255);">1</span>)*(*nrows)+<span style="color: rgb(0, 0, 255);">1</span>;<br />}<br /></pre></div>Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com0tag:blogger.com,1999:blog-1901659770749170850.post-78988674806021411342008-11-17T15:13:00.002-05:002008-11-18T13:07:21.465-05:00Call C from R and R from CSeveral years ago, while a research associate at the University of Chicago, I had the privilege of sitting in on a course taught by <a href="http://gsbportal.chicagogsb.edu/portal/server.pt/gateway/PTARGS_0_0_314_215_0_43/http%3B/gsbjob.chicagogsb.edu/Facultycourse08/FacultyDetail2.aspx?min_year=20084&max_year=20093&person_id=163822&fromSam=yes">Peter Rossi</a>: <a href="http://gsbportal.chicagogsb.edu/portal/server.pt/gateway/PTARGS_0_0_306_205_0_43/http%3B/gsbjob.chicagogsb.edu/Facultycourse08/CourseProfAllDetail.aspx?course_ns=37904&section_id=all&ac_year=2008">Bayesian Applications in Marketing and MicroEconometrics</a>. This course -- one I recommend to anyone at U Chicago who is interested in statistics -- was an incredibly clear treatment of Bayesian statistics, but the aspect I appreciated most was Peter's careful demonstration of Bayesian theory and methods using R.<br /><br />One feature of R that I had not made use of up until that point was the ability to call compiled C and Fortran functions from within R (this makes loop-heavy Metropolis-Hastings samplers much, much faster). It turns out that you can also include the R libraries in C source code so that R functions (e.g., random number generators) can be easily accessed. The R-Cran website has an excellent tutorial on how to develop R extensions (<a href="http://cran.r-project.org/doc/manuals/R-exts.html#System-and-foreign-language-interfaces">here</a>), but I wanted to share an example Peter used in class because it is extremely brief, and for 95% of what I do, this is all I need.<br /><br />As Peter writes, this is an incredibly inefficient way of simulating from the chisquare distribution, but it demonstrates the point. His more extensive writeup is located <a href="http://faculty.chicagogsb.edu/peter.rossi/teaching/37904/Adding%20Functions%20Written%20in%20C%20to%20R.pdf">here</a>.<br /><br />Save the following as testfun.c:<br /><div style="overflow: auto;"><pre><span style="color: rgb(0, 128, 0);">#include <R.h></span><br /><span style="color: rgb(0, 128, 0);">#include <Rmath.h></span><br /><span style="color: rgb(0, 128, 0);">#include <math.h></span><br /><span style="color: rgb(128, 128, 128);"><i>/* Function written by Peter Rossi from the University of Chicago GSB */</i></span><br /><span style="color: rgb(128, 128, 128);"><i>/*http://faculty.chicagogsb.edu/peter.rossi/teaching/37904/Adding%20Functions%20Written%20in%20C%20to%20R.pdf */</i></span><br /><br /><span style="color: rgb(128, 128, 128);"><i>/* include standard C math library and R internal function declarations */</i></span><br /><span style="color: rgb(128, 0, 0);">void</span> mychisq(<span style="color: rgb(128, 0, 0);">double</span> *vec, <span style="color: rgb(128, 0, 0);">double</span> *chisq, <span style="color: rgb(128, 0, 0);">int</span> *nu)<br /> <span style="color: rgb(128, 128, 128);"><i>/* void means return nothing */</i></span><br />{<br /> <span style="color: rgb(128, 0, 0);">int</span> i,iter; <span style="color: rgb(128, 128, 128);"><i>/* declare local vars */</i></span><br /> <span style="color: rgb(128, 128, 128);"><i>/* all statements end in; */</i></span><br /> GetRNGstate(); <span style="color: rgb(128, 128, 128);"><i>/* set random number seed */</i></span><br /> <b>for</b> (i=<span style="color: rgb(0, 0, 255);">0</span> ; i < *nu; ++i)<br /> <span style="color: rgb(128, 128, 128);"><i>/* loop over elements of vec */</i></span><br /> <span style="color: rgb(128, 128, 128);"><i>/*nu "dereferences" the pointer */</i></span><br /> { <span style="color: rgb(128, 128, 128);"><i>/* vectors start at 0 location!*/</i></span><br /> vec[i] = rnorm(<span style="color: rgb(128, 0, 128);">0.0</span>,<span style="color: rgb(128, 0, 128);">1.0</span>); <span style="color: rgb(128, 128, 128);"><i>/*use R function to draw normals */</i></span><br /> Rprintf(<span style="color: rgb(221, 0, 0);">"%ith normal draw= %lf </span><span style="color: rgb(255, 0, 255);">\n</span><span style="color: rgb(221, 0, 0);">"</span>,(i+<span style="color: rgb(0, 0, 255);">1</span>),vec[i]);<br /> <span style="color: rgb(128, 128, 128);"><i>/* print out results for "debugging" */</i></span><br /> }<br /> *chisq=<span style="color: rgb(128, 0, 128);">0.0</span>;<br /> iter=<span style="color: rgb(0, 0, 255);">0</span>;<br /> <b>while</b>(iter < *nu) <span style="color: rgb(128, 128, 128);"><i>/* "while" version of a loop */</i></span><br /> {<br /> <b>if</b>( iter == iter)<br /> {*chisq=*chisq + vec[iter]*vec[iter];}<br /> <span style="color: rgb(128, 128, 128);"><i>/* redundant if stmnt */</i></span><br /> iter=iter+<span style="color: rgb(0, 0, 255);">1</span>; <span style="color: rgb(128, 128, 128);"><i>/* note: can't use ** */</i></span><br /> <span style="color: rgb(128, 128, 128);"><i>/* if you want to be "cool" use iter += 1 */</i></span><br /> }<br /> PutRNGstate(); <span style="color: rgb(128, 128, 128);"><i>/* write back ran number seed */</i></span><br />}<br /><br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br /><br />To call this function in R, you first need to compile it. To do this you need all the standard compilers and libraries for your operating system. For Debian or Ubuntu, this should do it (if I missed a package, let me know in the comments):<br /><div style="overflow: auto;"><pre><br />$ sudo aptitude update<br />$ sudo aptitude install build-essential r-base-dev<br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br />Now, you should be able to compile the function:<br /><div style="overflow: auto;"><pre><br />$ R CMD SHLIB testfun.c<br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br />If all goes well, you should see the files testfun.o and testfun.so in the directory. To test the function we will source the following R script into R:<br /><div style="overflow: auto;"><pre><br />call_mychisq<span style="color: rgb(51, 99, 102);"><b><-</b></span><span style="color: rgb(0, 0, 191);">function</span>(nu){<br /><span style="color: rgb(176, 0, 0);"><b>##This function is just a wrapper for .C</b></span><br />vector<span style="color: rgb(51, 99, 102);"><b>=</b></span><b>double</b>(nu); chisq<span style="color: rgb(51, 99, 102);"><b>=</b></span><span style="color: rgb(0, 0, 255);">1</span><br /><b>.C</b>(<span style="color: rgb(221, 0, 0);">"mychisq"</span>,<b>as.double</b>(vector),<span style="color: rgb(128, 0, 0);">res=</span><b>as.double</b>(chisq),<b>as.integer</b>(nu))<span style="color: rgb(128, 63, 0);">$</span>res<br />}<br /><br /><span style="color: rgb(176, 0, 0);"><b>##Load the compiled code (you may need to include </b></span><br /><span style="color: rgb(176, 0, 0);"><b>## the explicit file path if it is not local<br />## NOTE: for Windows machines, you will want to load testfun.dll" </b></span><br /><b>dyn.load</b>(<span style="color: rgb(221, 0, 0);">"testfun.so"</span>)<br />result<span style="color: rgb(51, 99, 102);"><b><-</b></span><b>call_mychisq</b>(<span style="color: rgb(0, 0, 255);">10</span>)<br /><br /><span style="color: rgb(176, 0, 0);"><b>##If you re-compile testfun.c, you must unload it</b></span><br /><span style="color: rgb(176, 0, 0);"><b>## and then re-load it:</b></span><br /><span style="color: rgb(176, 0, 0);"><b>## dyn.unload("testfun.so")</b></span><br /><span style="color: rgb(176, 0, 0);"><b>## dyn.load("testfun.so")</b></span><br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br />And get the following output:<br /><div style="overflow: auto;"><pre><br />> dyn.load("testfun.so")<br />> result<-call_mychisq(10) 1th normal draw= -1.031170 2th normal draw= -1.214103 3th normal draw= 0.002335 4th normal draw= 0.296146 5th normal draw= -0.908862 6th normal draw= -1.567820 7th normal draw= -0.079227 8th normal draw= -1.404414 9th normal draw= 0.616567 10th normal draw= -0.007855 > result<br />[1] 8.268028<br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div>Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com3tag:blogger.com,1999:blog-1901659770749170850.post-84763660030425316212008-11-05T14:47:00.003-05:002009-08-21T10:40:47.355-04:00Using subversion to manage codeI have finally come to terms with the fact that I need some kind of version control for the projects I am working on and the best bet these days is Subversion (svn). I have been using svn for some time now via a GUI client (Linux: kdesvn, Windows: tortoisesvn); however, it turns out that working with svn from the command line is pretty easy and far more versatile. For a complete treatment of this subject, check out the online documentation <a href="http://svnbook.red-bean.com/en/1.1/index.html">here</a> and an excellent cheat sheet <a href="http://www.abbeyworkshop.com/howto/misc/svn01/">here</a>. What follows is a quick primer on the very basics of setting up and managing your svn repositories on a local machine or server.<br /><br />For ease of use I will describe the creation of a single repository for a single project. This means a little more overhead; however, it makes the repository more portable and flexible in the long run.<br /><br />1) First we set up the repository structure in a temporary folder (either on the server or locally):<br /><div style="overflow: auto;"><pre>$ mkdir ~/tmp<br />$ mkdir ~/tmp/project1<br />$ mkdir ~/tmp/project1/trunk<br />$ mkdir ~/tmp/project1/branches<br />$ mkdir ~/tmp/project1/tags</pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br />2) Now, make a folder to hold your repositories and create an empty repository for your project.<br /><div style="overflow: auto;"><pre>$ mkdir ~/svnroot<br />$ svnadmin create ~/svnroot/project1</pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br />3) Import the folder structure into the empty repository. After this import, the folders in tmp can be removed -- they are only there to make the creation of the folder structure easier.<br /><div style="overflow: auto;"><pre>$ svn import ~/tmp/project1 file:///home/myusername/svnroot/project1 -m "Initial import"</pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br />4) Finally, make your current project folder a "working copy" of the repository. Checkout the trunk (or head) of the repository to the folder where project1 currently resides (in this example, the existing project files are located at ~/working/project1).<br />If you created your repository on a local folder:<br /><div style="overflow: auto;"><pre>$ svn checkout file:///home/myusername/svnroot/project1/trunk ~/working/project1</pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br />Alternatively, if you created your repository on a remote server:<br /><div style="overflow: auto;"><pre>$ svn checkout svn+ssh://remote.server.name/home/myusername/svnroot/project1/trunk ~/working/project1</pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br />Because the repository is empty at this stage, all the above commands do is create a .svn folder in the ~/working/project1 directory. The following command will show that there are folders and files in the project directory that are not currently part of the repository:<br /><div style="overflow: auto;"><pre>$ svn st<br />? somefolder<br />? someotherfolder<br />? somefile.txt</pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br />Now you need to add all of the files and folders in this directory to the repository. This is easily accomplished using a bit of awk code (<a href="http://codesnippets.joyent.com/posts/show/935#comments">modified from a post here</a>):<br /><br /><div style="overflow: auto;"><pre>$ svn status | grep "^\?" | awk -F " " '{print $2}' | tr "\n" "\0" | xargs -0 svn add<br />$ svn st<br />A somefolder<br />A somefolder/file1.txt<br />A somefolder/file2.txt<br />A someotherfolder<br />A someotherfolder/file3.txt<br />A somefile.txt</pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br />Now you just need to commit these changes and your working directory is up to date:<br /><div style="overflow: auto;"><pre>$ svn commit -m "Adding original files to repository."</pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br /><br />5) When you are ready to commit new changes to the repository, make sure that all new files/folders are added and all deleted files/folders are removed:<br /><div style="overflow: auto;"><pre>$ svn status | grep "^\?" | awk -F " " '{print $2}' | tr "\n" "\0" | xargs -0 svn add<br />$ svn status | grep "^\!" | awk -F " " '{print $2}' | tr "\n" "\0" | xargs -0 svn remove<br />$ svn commit -m "Some comment to remind you why you are committing changes..."</pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br /><br />6) Finally, if for some reason you want to remove a working directory from versioning (i.e., get rid of the .svn folders that are placed in every folder subfolder), use the following:<br /><div style="overflow: auto;"><pre>$ cd ~/working/project1<br />$ rm -rf 'find . -name .svn'</pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br /><br />Update: If your svn repository changes to a new server name, use the following syntax to update your working directory:<br /><div style="overflow: auto;"><pre>$ cd ~/working/project1<br />$svn info<br />[the current URL and other info are printed to the screen]<br />$svn switch --relocate svn+ssh://OLD.URL/path/to/svnrepo svn+ssh://NEW.URL/path/to/svnrepo .<br />$svn commit -m "new update"<br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br />Note: If you mis-type the old URL, "svn switch" will fail silently. So make sure to check that it has updated by using the "svn info" command.Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com0tag:blogger.com,1999:blog-1901659770749170850.post-35934993817725833652008-09-20T15:07:00.000-04:002008-09-21T13:56:44.289-04:00Installing Ubuntu Linux on existing LVM partition.I wanted to install Ubuntu on a computer that already housed Fedora Linux. It seems that Fedora defaults to installing /root and /swap on an <a href="http://en.wikipedia.org/wiki/Logical_volume_management">LVM partition</a> -- something that Ubuntu does not recognize out of the box. After poking about online, here is the process I followed (compiled from <a href="http://www.linux-sxs.org/storage/fedora2ubuntu.html">here</a>, <a href="http://www.faqs.org/docs/Linux-HOWTO/LVM-HOWTO.html#AEN638">here</a>, and <a href="http://www.debuntu.org/how-to-install-ubuntu-over-lvm-filesystem">here</a>.):<br /><br /><div style="overflow: auto; height: 400px;"><br /><span style="color: rgb(255,0,0)">Warning: As always, any time you are messing about with drive partitions, you risk something tragic happening to your data. Backup critical bits before continuing.</span><br /><br />First deal with reducing the existing logical volume:<br /><blockquote><br /># Boot Ubuntu from live-CD.<br /><br /># Install lvm2:<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">$ sudo apt-get install lvm2<br /></span><br /># Load the necessary module(s):<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">$ sudo modprobe dm-mod</span><br /><br /># Scan your system for LVM volumes and identify in the output the volume group name that has your Fedora volume (mine proved to be VolGroup00):<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">$ sudo vgscan</span><br /><br /># Activate the volume:<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">$ sudo vgchange -ay VolGroup00</span><br /><br /># Find the logical volume that has your Fedora root filesystem (mine proved to be LogVol00):<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">$ sudo lvs</span><br /><br /># Determine how much space is used on the volume:<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">$ sudo vgdisplay</span><br /><br /># Check the filesystem<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">sudo e2fsck -f /dev/VolGroup00/LogVol00</span><br /><br /># Resize the filesystem to desired size (CRITICAL before proceeding to next step!)<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">$ resize2fs -f /dev/VolGroup00/LogVol00 16G</span><br /><br /># Reduce the size of the logical volume<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">$ lvreduce -L10G /dev/VolGroup00/LogVol00</span></blockquote><br /><br />Create new logical volumes for /root and /home (in my case /swap already existed):<br /><blockquote><br /># First make sure how much free space remains in the volume group<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">$ sudo vgdisplay</span><br /><br /># Create \UbuntuRoot and \home<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">$ sudo lvcreate -n UbuntuRoot -L10G VolGroup00<br />$ sudo lvcreate -n home -L18.28G VolGroup00</span><br /><br /># Now format the new logical volumes:<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">$ sudo mkfs -j /dev/VolGroup00/UbuntuRoot -L root<br />$ sudo mkfs -j /dev/VolGroup00/home -L home</span><br /></blockquote><br /><br />Then run the installer, choosing the appropriate volumes. Note that I had already created a hard partition to serve as /boot --- this is the only way a lvm install will work.<br /><br />Finally, make sure that the newly installed system has lvm2 included (otherwise it will not start). The following is adapted from <a href="http://www.debuntu.org/how-to-install-ubuntu-over-lvm-filesystem">here</a>.<br /><br /><blockquote><br /># Create a mountpoint for the newly installed system:<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">$ sudo mkdir /target<br />$ sudo mkdir /target</span><br /><br /># Mount the newly installed system:<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">$ sudo mount /dev/VolGroup00/UbuntuRoot /target<br />$ sudo mount /dev/sda6 /target/boot</span><br /><br /># Activate those directories and install lvm2<br /><span style="background-color: rgb(200, 200, 200); color: rgb(0, 0, 0);">$ sudo chroot /target<br />$ sudo apt-get update<br />$ sudo apt-get install lvm2</span><br /></blockquote><br /><br />Now you can restart into the new system!<br /><br /><br /><br /></div>Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com0tag:blogger.com,1999:blog-1901659770749170850.post-33525662518138718302008-05-06T10:21:00.000-04:002008-05-06T10:53:06.811-04:00Error captureIn a recent post to r-sig-ecology, Mike Colvin suggested the following to capture errors within a loop:<br /><div style="overflow: auto;"><pre><br /><span style="color: #0000bf;">for</span><span style="color: #000000;"> </span><span style="font-weight: bold;color: #000000;">(</span><span style="color: #000000;">i </span><span style="color: #0000bf;">in</span><span style="color: #000000;"> </span><span style="color: #0000ff;">1</span><span style="color: #803f00;">:</span><span style="color: #0000ff;">1000</span><span style="font-weight: bold;color: #000000;">)</span><span style="color: #000000;">{</span><br /><span style="color: #000000;">fit</span><span style="font-weight: bold;color: #336366;"><-</span><span style="font-weight: bold;color: #000000;">try(lm(</span><span style="color: #000000;">y</span><span style="color: #803f00;">~</span><span style="color: #000000;">x,dataset</span><span style="font-weight: bold;color: #000000;">))</span><br /><span style="color: #000000;">results</span><span style="font-weight: bold;color: #336366;"><-</span><span style="color: #000000;"> </span><span style="font-weight: bold;color: #000000;">ifelse(class(</span><span style="color: #000000;">fit</span><span style="font-weight: bold;color: #000000;">)</span><span style="color: #803f00;">==</span><span style="color: #dd0000;">"try-error"</span><span style="color: #000000;">, </span><span style="color: #008000;">NA</span><span style="color: #000000;">, fit</span><span style="color: #803f00;">$</span><span style="color: #000000;">coefficients</span><span style="font-weight: bold;color: #000000;">)</span><br /><span style="color: #000000;">}</span><br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div>Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com0tag:blogger.com,1999:blog-1901659770749170850.post-38385317010004309752008-03-18T15:04:00.003-04:002011-01-21T15:38:55.527-05:00Plotting contours<a onblur="try {parent.deselectBloggerImageGracefully();} catch(e) {}" href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEikJ3OFvB8g9p7yPBFavHShbuLGJkYhW5Kifoa2mrraD2j5GMa-jahqHQD79fw-0EhssA_SGwdFGWOvhc4Y2JHO_yh1LHV1V4videGkAyFf9dXh1UJyRZa7I60HWxq_MPQcwcAQfBNFxWw4/s1600-h/densplot.png"><img style="margin: 0px auto 10px; display: block; text-align: center; cursor: pointer;" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEikJ3OFvB8g9p7yPBFavHShbuLGJkYhW5Kifoa2mrraD2j5GMa-jahqHQD79fw-0EhssA_SGwdFGWOvhc4Y2JHO_yh1LHV1V4videGkAyFf9dXh1UJyRZa7I60HWxq_MPQcwcAQfBNFxWw4/s320/densplot.png" alt="" id="BLOGGER_PHOTO_ID_5179188337696722914" border="0" /></a><br />Plenty of packages allow you to plot contours of a "z" value; however, I wanted to be able to plot a specific density contour of a sample from a bivariate distribution over a plot that was a function of the x and y parameters. The example only plots the density of x and y; however, if you set plot.dens=FALSE, the contours will be added to the active display device (i.e., over an image(x,y,function(x,y){...})).<br /><br />This function also lets you specify the line widths and types for each of the contours.<br /><div style="overflow: auto;"><pre><span style="font-weight: bold; color: rgb(176, 0, 0);">###########################################</span><span style="font-weight: bold; color: rgb(176, 0, 0);"><br />## drawcontour.R</span><br /><span style="font-weight: bold; color: rgb(176, 0, 0);">## Written by J.D. Forester, 17 March 2008</span><br /><span style="font-weight: bold; color: rgb(176, 0, 0);">###########################################</span><br /><pre><br /><span style="font-weight: bold; color: rgb(176, 0, 0);">##This function draws an approximate density contour based on</span><br /><span style="font-weight: bold; color: rgb(176, 0, 0);">##empirical, bivariate data.</span><br /><br /><br /><span style="font-weight: bold; color: rgb(176, 0, 0);">##change testit to FALSE if sourcing the file</span><br /><span style="color: rgb(0, 0, 0);">testit</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="color: rgb(0, 128, 0);">TRUE</span><br /><br /><br /><span style="color: rgb(0, 0, 0);">draw.contour</span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="color: rgb(0, 0, 191);">function</span><span style="font-weight: bold; color: rgb(0, 0, 0);">(</span><span style="color: rgb(0, 0, 0);">a,</span><span style="color: rgb(128, 0, 0);">alpha=</span><span style="color: rgb(128, 0, 128);">0.95</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(128, 0, 0);">plot.dens=</span><span style="color: rgb(0, 128, 0);">FALSE</span><span style="color: rgb(0, 0, 0);">, </span><span style="color: rgb(128, 0, 0);">line.width=</span><span style="color: rgb(0, 0, 255);">2</span><span style="color: rgb(0, 0, 0);">, </span><span style="color: rgb(128, 0, 0);">line.type=</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">, </span><span style="color: rgb(128, 0, 0);">limits=</span><span style="color: rgb(0, 128, 0);">NULL</span><span style="color: rgb(0, 0, 0);">, </span><span style="color: rgb(128, 0, 0);">density.res=</span><span style="color: rgb(0, 0, 255);">300</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(128, 0, 0);">spline.smooth=</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">,...</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">{</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##a is a list or matrix of x and y coordinates (e.g., a=list("x"=rnorm(100),"y"=rnorm(100)))</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">## if a is a list or dataframe, the components must be labeled "x" and "y"</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">## if a is a matrix, the first column is assumed to be x, the second y</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##alpha is the contour level desired</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##if plot.dens==TRUE, then the joint density of x and y are plotted,</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">## otherwise the contour is added to the current plot.</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##density.res controls the resolution of the density plot</span><br /><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##A key assumption of this function is that very little probability mass lies outside the limits of</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">## the x and y values in "a". This is likely reasonable if the number of observations in a is large.</span><br /><span style="color: rgb(0, 0, 0);"> </span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">require(</span><span style="color: rgb(0, 0, 0);">MASS</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">require(</span><span style="color: rgb(0, 0, 0);">ks</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 191);">if</span><span style="font-weight: bold; color: rgb(0, 0, 0);">(length(</span><span style="color: rgb(0, 0, 0);">line.width</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(128, 63, 0);">!=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">length(</span><span style="color: rgb(0, 0, 0);">alpha</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><span style="color: rgb(0, 0, 0);">{</span><br /><span style="color: rgb(0, 0, 0);"> line.width </span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">rep(</span><span style="color: rgb(0, 0, 0);">line.width[</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">],</span><span style="font-weight: bold; color: rgb(0, 0, 0);">length(</span><span style="color: rgb(0, 0, 0);">alpha</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><br /><span style="color: rgb(0, 0, 0);"> }</span><br /><br /><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 191);">if</span><span style="font-weight: bold; color: rgb(0, 0, 0);">(length(</span><span style="color: rgb(0, 0, 0);">line.type</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(128, 63, 0);">!=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">length(</span><span style="color: rgb(0, 0, 0);">alpha</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><span style="color: rgb(0, 0, 0);">{</span><br /><span style="color: rgb(0, 0, 0);"> line.type </span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">rep(</span><span style="color: rgb(0, 0, 0);">line.type[</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">],</span><span style="font-weight: bold; color: rgb(0, 0, 0);">length(</span><span style="color: rgb(0, 0, 0);">alpha</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><br /><span style="color: rgb(0, 0, 0);"> }</span><br /><span style="color: rgb(0, 0, 0);"> </span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 191);">if</span><span style="font-weight: bold; color: rgb(0, 0, 0);">(is.matrix(</span><span style="color: rgb(0, 0, 0);">a</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><span style="color: rgb(0, 0, 0);">{</span><br /><span style="color: rgb(0, 0, 0);"> a</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">list(</span><span style="color: rgb(221, 0, 0);">"x"</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="color: rgb(0, 0, 0);">a[,</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">],</span><span style="color: rgb(221, 0, 0);">"y"</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="color: rgb(0, 0, 0);">a[,</span><span style="color: rgb(0, 0, 255);">2</span><span style="color: rgb(0, 0, 0);">]</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><span style="color: rgb(0, 0, 0);"> }</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##generate approximate density values</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 191);">if</span><span style="font-weight: bold; color: rgb(0, 0, 0);">(is.null(</span><span style="color: rgb(0, 0, 0);">limits</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><span style="color: rgb(0, 0, 0);">{</span><br /><span style="color: rgb(0, 0, 0);"> limits</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">c(range(</span><span style="color: rgb(0, 0, 0);">a</span><span style="color: rgb(128, 63, 0);">$</span><span style="color: rgb(0, 0, 0);">x</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">,</span><span style="font-weight: bold; color: rgb(0, 0, 0);">range(</span><span style="color: rgb(0, 0, 0);">a</span><span style="color: rgb(128, 63, 0);">$</span><span style="color: rgb(0, 0, 0);">y</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><br /><span style="color: rgb(0, 0, 0);"> }</span><br /><span style="color: rgb(0, 0, 0);"> f1</span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="font-weight: bold; color: rgb(0, 0, 0);">kde2d(</span><span style="color: rgb(0, 0, 0);">a</span><span style="color: rgb(128, 63, 0);">$</span><span style="color: rgb(0, 0, 0);">x,a</span><span style="color: rgb(128, 63, 0);">$</span><span style="color: rgb(0, 0, 0);">y,</span><span style="color: rgb(128, 0, 0);">n=</span><span style="color: rgb(0, 0, 0);">density.res,</span><span style="color: rgb(128, 0, 0);">lims=</span><span style="color: rgb(0, 0, 0);">limits</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##plot empirical density</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 191);">if</span><span style="font-weight: bold; color: rgb(0, 0, 0);">(</span><span style="color: rgb(0, 0, 0);">plot.dens</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">image(</span><span style="color: rgb(0, 0, 0);">f1,...</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><span style="color: rgb(0, 0, 0);"> </span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 191);">if</span><span style="font-weight: bold; color: rgb(0, 0, 0);">(is.null(dev.list()))</span><span style="color: rgb(0, 0, 0);">{</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##ensure that there is a window in which to draw the contour</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">plot(</span><span style="color: rgb(0, 0, 0);">a,</span><span style="color: rgb(128, 0, 0);">type=</span><span style="color: rgb(221, 0, 0);">"n"</span><span style="color: rgb(221, 0, 0);">,</span><span style="color: rgb(128, 0, 0);">xlim=</span><span style="color: rgb(0, 0, 0);">limits[1:2]</span><span style="color: rgb(221, 0, 0);">,</span><span style="color: rgb(128, 0, 0);">ylim=</span><span style="color: rgb(0, 0, 0);">limits[3:4]</span><span style="font-weight: bold; color: rgb(0, 0, 0);">,...)</span><br /><span style="color: rgb(0, 0, 0);"> }</span><br /><span style="color: rgb(0, 0, 0);"> </span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##estimate critical contour value</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">## assume that density outside of plot is very small</span><br /><br /><span style="color: rgb(0, 0, 0);"> zdens </span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">rev(sort(</span><span style="color: rgb(0, 0, 0);">f1</span><span style="color: rgb(128, 63, 0);">$</span><span style="color: rgb(0, 0, 0);">z</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><br /><span style="color: rgb(0, 0, 0);"> Czdens </span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">cumsum(</span><span style="color: rgb(0, 0, 0);">zdens</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><span style="color: rgb(0, 0, 0);"> Czdens </span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">(</span><span style="color: rgb(0, 0, 0);">Czdens</span><span style="color: rgb(128, 63, 0);">/</span><span style="color: rgb(0, 0, 0);">Czdens[</span><span style="font-weight: bold; color: rgb(0, 0, 0);">length(</span><span style="color: rgb(0, 0, 0);">zdens</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">]</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 191);">for</span><span style="font-weight: bold; color: rgb(0, 0, 0);">(</span><span style="color: rgb(0, 0, 0);">cont.level </span><span style="color: rgb(0, 0, 191);">in</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(128, 63, 0);">:</span><span style="font-weight: bold; color: rgb(0, 0, 0);">length(</span><span style="color: rgb(0, 0, 0);">alpha</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><span style="color: rgb(0, 0, 0);">{</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##This loop allows for multiple contour levels</span><br /><span style="color: rgb(0, 0, 0);"> crit.val </span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="color: rgb(0, 0, 0);"> zdens[</span><span style="font-weight: bold; color: rgb(0, 0, 0);">max(which(</span><span style="color: rgb(0, 0, 0);">Czdens</span><span style="color: rgb(128, 63, 0);"><=</span><span style="color: rgb(0, 0, 0);">alpha[cont.level]</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><span style="color: rgb(0, 0, 0);">]</span><br /><span style="color: rgb(0, 0, 0);"> </span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##determine coordinates of critical contour</span><br /><span style="color: rgb(0, 0, 0);"> b.full</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">contourLines(</span><span style="color: rgb(0, 0, 0);">f1,</span><span style="color: rgb(128, 0, 0);">levels=</span><span style="color: rgb(0, 0, 0);">crit.val</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 191);">for</span><span style="font-weight: bold; color: rgb(0, 0, 0);">(</span><span style="color: rgb(0, 0, 0);">c </span><span style="color: rgb(0, 0, 191);">in</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(128, 63, 0);">:</span><span style="font-weight: bold; color: rgb(0, 0, 0);">length(</span><span style="color: rgb(0, 0, 0);">b.full</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><span style="color: rgb(0, 0, 0);">{</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##This loop is used in case the density is multimodal or if the desired contour</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">## extends outside the plotting region</span><br /><span style="color: rgb(0, 0, 0);"> b</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">list(</span><span style="color: rgb(221, 0, 0);">"x"</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">as.vector(unlist(</span><span style="color: rgb(0, 0, 0);">b.full[[c]][</span><span style="color: rgb(0, 0, 255);">2</span><span style="color: rgb(0, 0, 0);">]</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(221, 0, 0);">"y"</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">as.vector(unlist(</span><span style="color: rgb(0, 0, 0);">b.full[[c]][</span><span style="color: rgb(0, 0, 255);">3</span><span style="color: rgb(0, 0, 0);">]</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)))</span><br /><span style="color: rgb(0, 0, 0);"> </span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##plot desired contour</span><br /><span style="color: rgb(0, 0, 0);"> line.dat</span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="font-weight: bold; color: rgb(0, 0, 0);">xspline(</span><span style="color: rgb(0, 0, 0);">b,</span><span style="color: rgb(128, 0, 0);">shape=</span><span style="color: rgb(0, 0, 0);">spline.smooth,</span><span style="color: rgb(128, 0, 0);">open=</span><span style="color: rgb(0, 128, 0);">TRUE</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(128, 0, 0);">draw=</span><span style="color: rgb(0, 128, 0);">FALSE</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">lines(</span><span style="color: rgb(0, 0, 0);">line.dat,</span><span style="color: rgb(128, 0, 0);">lty=</span><span style="color: rgb(0, 0, 0);">line.type[cont.level],</span><span style="color: rgb(128, 0, 0);">lwd=</span><span style="color: rgb(0, 0, 0);">line.width[cont.level]</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><span style="color: rgb(0, 0, 0);"> }</span><br /><span style="color: rgb(0, 0, 0);"> }</span><br /><span style="color: rgb(0, 0, 0);">}</span><br /><br /><span style="font-weight: bold; color: rgb(176, 0, 0);">##############################</span><br /><span style="font-weight: bold; color: rgb(176, 0, 0);">##Test the function</span><br /><span style="font-weight: bold; color: rgb(176, 0, 0);">##############################</span><br /><br /><span style="font-weight: bold; color: rgb(176, 0, 0);">##generate data</span><br /><span style="color: rgb(0, 0, 191);">if</span><span style="font-weight: bold; color: rgb(0, 0, 0);">(</span><span style="color: rgb(0, 0, 0);">testit</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">{</span><br /><span style="color: rgb(0, 0, 0);"> n</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="color: rgb(0, 0, 255);">10000</span><br /><span style="color: rgb(0, 0, 0);"> a</span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="font-weight: bold; color: rgb(0, 0, 0);">list(</span><span style="color: rgb(221, 0, 0);">"x"</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">rnorm(</span><span style="color: rgb(0, 0, 0);">n,</span><span style="color: rgb(0, 0, 255);">400</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(0, 0, 255);">100</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(221, 0, 0);">"y"</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">rweibull(</span><span style="color: rgb(0, 0, 0);">n,</span><span style="color: rgb(0, 0, 255);">2</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(0, 0, 255);">100</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">draw.contour(</span><span style="color: rgb(128, 0, 0);">a=</span><span style="color: rgb(0, 0, 0);">a,</span><span style="color: rgb(128, 0, 0);">alpha=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">c(</span><span style="color: rgb(128, 0, 128);">0.95</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(128, 0, 128);">0.5</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(128, 0, 128);">0.05</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(128, 0, 0);">line.width=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">c(</span><span style="color: rgb(0, 0, 255);">2</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(0, 0, 255);">2</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(128, 0, 0);">line.type=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">c(</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(0, 0, 255);">2</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(0, 0, 255);">1</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(128, 0, 0);">plot.dens=</span><span style="color: rgb(0, 128, 0);">TRUE</span><span style="font-weight: bold; color: rgb(0, 0, 0);">,</span><span style="color: rgb(128, 0, 0);"> xlab=</span>"X"<span style="font-weight: bold; color: rgb(0, 0, 0);">,<span style="font-family:monospace;"> </span></span><span style="color: rgb(128, 0, 0);">ylab=</span>"Y"<span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><span style="color: rgb(0, 0, 0);">}</span><br /></pre><br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div>Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com3tag:blogger.com,1999:blog-1901659770749170850.post-41956409452801642862008-02-04T11:44:00.004-05:002008-12-01T11:03:49.286-05:00Drop unused factor levelsWhen creating a subset of a dataframe, I often exclude rows based on the level of a factor. However, the "levels" of the factor remain intact. This is the intended behavior of R, but it can cause problems in some cases. I finally discovered how to clean up levels in <a href="http://www.nabble.com/Bug-in-levels%28%29-function--td15131875.html#a15135262">this post to R-Help</a>. Here is an example:<br /><div style="overflow: auto;"><pre><span style="color: rgb(0, 0, 0);">> a </span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">factor(</span><span style="color: rgb(0, 0, 0);">letters</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><span style="color: rgb(0, 0, 0);">> a</span><br /><span style="color: rgb(0, 0, 0);"> [</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">] a b c d e f g h i j k l m n o p q r s t u v w x y z</span><br /><span style="color: rgb(0, 0, 0);">Levels</span><span style="color: rgb(128, 63, 0);">:</span><span style="color: rgb(0, 0, 0);"> a b c d e f g h i j k l m n o p q r s t u v w x y z</span><br /><br /><span style="font-weight: bold; color: rgb(176, 0, 0);">## Now, even though b only includes five letters, </span><br /><span style="font-weight: bold; color: rgb(176, 0, 0);">## all 26 are listed in the levels</span><br /><span style="color: rgb(0, 0, 0);">> b </span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="color: rgb(0, 0, 0);"> a[</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(128, 63, 0);">:</span><span style="color: rgb(0, 0, 255);">5</span><span style="color: rgb(0, 0, 0);">]</span><br /><span style="color: rgb(0, 0, 0);">> b</span><br /><span style="color: rgb(0, 0, 0);">[</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">] a b c d e</span><br /><span style="color: rgb(0, 0, 0);">Levels</span><span style="color: rgb(128, 63, 0);">:</span><span style="color: rgb(0, 0, 0);"> a b c d e f g h i j k l m n o p q r s t u v w x y z</span><br /><br /><span style="font-weight: bold; color: rgb(176, 0, 0);">## This behavior can be changed using the following syntax:</span><br /><span style="color: rgb(0, 0, 0);">> b </span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="color: rgb(0, 0, 0);"> a[</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(128, 63, 0);">:</span><span style="color: rgb(0, 0, 255);">5</span><span style="color: rgb(0, 0, 0);">,drop </span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 128, 0);">TRUE</span><span style="color: rgb(0, 0, 0);">]</span><br /><span style="color: rgb(0, 0, 0);">> b</span><br /><span style="color: rgb(0, 0, 0);">[</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">] a b c d e</span><br /><span style="color: rgb(0, 0, 0);">Levels</span><span style="color: rgb(128, 63, 0);">:</span><span style="color: rgb(0, 0, 0);"> a b c d e</span><br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br />Another way to deal with this is to use the dropUnusedLevels() command in the Hmisc library. The only issue here is that behavior is changed globally which may have undesired consequences (see the post listed above).<br /><br /><span style="font-weight: bold;">****UPDATE****</span><br />As Jeff Hollister mentions in the comments, there is another way to do this:<br /><div style="overflow:auto"><pre><br />a<-factor(letters)<br />b<-factor(a[1:5])<br /></pre><span style="font-size:85%;"><span style="font-family:courier<br />new;"></span></span></div><br /><br />Yet another way, if you are working with data frames that by default convert characters into factors, was suggested on r-sig-ecology by Hadley Wickham:<br /><div style="overflow:auto"><pre><br />options(stringsAsFactors = FALSE)<br />a <-data.frame("alpha"=letters)<br />b<-a[1:5]<br /></pre><span style="font-size:85%;"><span style="font-family:courier<br />new;"></span></span></div>Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com4tag:blogger.com,1999:blog-1901659770749170850.post-22806866795712597682007-11-29T14:11:00.001-05:002009-08-31T09:54:34.100-04:00Convert factors to numbersIf you have a vector of factors it is easy to get the factor level; however, I always forget how to extract the factor value. I ran into the answer <a href="https://stat.ethz.ch/pipermail/r-help/2000-April/006423.html">here</a>.<br /><div id="thumbFrame" style="overflow: auto;"><pre><br /><span style="color: rgb(0, 0, 0);">> x</span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="font-weight: bold; color: rgb(0, 0, 0);">factor(c(round(rnorm(</span><span style="color: rgb(0, 0, 255);">10</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(0, 0, 255);">2</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(221, 0, 0);">"A"</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(221, 0, 0);">"B"</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(0, 128, 0);">NA</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><br /><span style="color: rgb(0, 0, 0);">> x</span><br /><span style="color: rgb(0, 0, 0);"> [</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">] </span><span style="color: rgb(128, 0, 128);">1.61</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">1.12</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">1.26</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">0.09</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(128, 0, 128);">0.13</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">0.16</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(128, 0, 128);">0.03</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(128, 0, 128);">0.1</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">0.09</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(128, 0, 128);">0.47</span><span style="color: rgb(0, 0, 0);"> </span><br /><span style="color: rgb(0, 0, 0);"> [</span><span style="color: rgb(0, 0, 255);">11</span><span style="color: rgb(0, 0, 0);">] A B </span><br /><span style="color: rgb(0, 0, 0);">Levels</span><span style="color: rgb(128, 63, 0);">:</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(128, 0, 128);">0.03</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">0.09</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(128, 0, 128);">0.1</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(128, 0, 128);">0.13</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">0.16</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(128, 0, 128);">0.47</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">1.12</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">1.26</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">1.61</span><span style="color: rgb(0, 0, 0);"> A B</span><br /><span style="color: rgb(0, 0, 0);">> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">as.numeric(</span><span style="color: rgb(0, 0, 0);">x</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><span style="color: rgb(0, 0, 0);"> [</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">] </span><span style="color: rgb(0, 0, 255);">9</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 255);">7</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 255);">8</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 255);">2</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 255);">4</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 255);">5</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 255);">3</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 255);">2</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 255);">6</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 255);">10</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 255);">11</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 128, 0);">NA</span><br /><span style="color: rgb(0, 0, 0);">> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">as.numeric(levels(</span><span style="color: rgb(0, 0, 0);">x</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">[x]</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><span style="color: rgb(0, 0, 0);"> [</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">] </span><span style="color: rgb(128, 0, 128);">1.61</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">1.12</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">1.26</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">0.09</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(128, 0, 128);">0.13</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">0.16</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(128, 0, 128);">0.03</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(128, 0, 128);">0.10</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 0, 128);">0.09</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(128, 0, 128);">0.47</span><span style="color: rgb(0, 0, 0);"> </span><br /><span style="color: rgb(0, 0, 0);"> [</span><span style="color: rgb(0, 0, 255);">11</span><span style="color: rgb(0, 0, 0);">] </span><span style="color: rgb(0, 128, 0);">NA</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 128, 0);">NA</span><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 128, 0);">NA</span><br /><span style="color: rgb(0, 0, 0);">Warning message</span><span style="color: rgb(128, 63, 0);">:</span><br /><span style="color: rgb(0, 0, 0);">NAs introduced by coercion </span></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span></div><br /><br /><br /><span style="font-weight: bold;">**UPDATE**</span><br />And another method mentioned in the comments (that I like better):<br /><br />> as.numeric(as.character(x))Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com1tag:blogger.com,1999:blog-1901659770749170850.post-5745420411002059642007-11-15T12:36:00.000-05:002008-09-21T12:06:04.398-04:00Preparing plots for publicationThe plotting capabilities of R are excellent; however, when I am preparing a figure for publication, I often need to combine multiple plots or add objects (e.g., arrows or text) to an existing plot. While this can be accomplished in R, my patience for tweaking layout parameters tends to run out quickly. I searched around and found a nice solution using open source software (an example with Matlab plots is described <a href="http://www.reimeika.ca/marco/matlab/coolplots.html">here</a>).<br /><br /><div style="overflow:auto"><br />1) Create the desired figure in R and export as an eps file:<br /><br /><br /><pre><br />library(ggplot2)<br />postscript("myfig.eps",horizontal=FALSE, paper="special",height=10,width=10)<br />qplot(x,y,data=mydat)<br />dev.off()<br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span><br /><br />2) Install the necessary software on your computer, and then convert the image file to one more easily edited.<br /><br /><pre><br />#install software<br />sudo aptitude install pstoedit tgif<br /><br />#clean up the eps file so it is more easily read<br />eps2eps myfig.eps myfig2.eps<br /><br />#convert the eps file to a tgif object<br />pstoedit -f tgif myfig2.eps myfig2.obj<br /></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span><br /><br />3) Edit the new file myfig2.obj in tgif.<br /><br /></div>Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com0tag:blogger.com,1999:blog-1901659770749170850.post-45014167885576933642007-10-18T00:27:00.001-04:002011-07-29T13:30:05.929-04:00Approximate sunrise and sunset timesThis function is not perfect, but it does a reasonable job estimating sunrise and sunset times for my field site. If more accurate data are required, try <a href="http://aa.usno.navy.mil/data/docs/RS_OneYear.php">here</a>. Note, the command to calculate day of year is: strptime(x, "%m/%d/%Y")$yday+1<br />
<div style="overflow: auto;"><pre><pre><code><span style="color: black;">suncalc</span><span style="color: #336366; font-weight: bold;"><-</span><span style="color: #0000bf;">function</span><span style="color: black; font-weight: bold;">(</span><span style="color: black;">d,</span><span style="color: maroon;">Lat=</span><span style="color: purple;">48.1442</span><span style="color: black;">,</span><span style="color: maroon;">Long=</span><span style="color: #803f00;">-</span><span style="color: purple;">122.7551</span><span style="color: black; font-weight: bold;">)</span><span style="color: black;">{</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## d is the day of year</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## Lat is latitude in decimal degrees</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## Long is longitude in decimal degrees (negative == West)</span>
<span style="color: black;"> </span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">##This method is copied from:</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">##Teets, D.A. 2003. Predicting sunrise and sunset times.</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## The College Mathematics Journal 34(4):317-321.</span>
<span style="color: black;"> </span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## At the default location the estimates of sunrise and sunset are within</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## seven minutes of the correct times (http://aa.usno.navy.mil/data/docs/RS_OneYear.php)</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## with a mean of 2.4 minutes error.</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## Function to convert degrees to radians</span>
<span style="color: black;"> rad</span><span style="color: #336366; font-weight: bold;"><-</span><span style="color: #0000bf;">function</span><span style="color: black; font-weight: bold;">(</span><span style="color: black;">x</span><span style="color: black; font-weight: bold;">)</span><span style="color: black;">pi</span><span style="color: #803f00;">*</span><span style="color: black;">x</span><span style="color: #803f00;">/</span><span style="color: blue;">180</span>
<span style="color: black;"> </span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">##Radius of the earth (km)</span>
<span style="color: black;"> R</span><span style="color: #336366; font-weight: bold;">=</span><span style="color: blue;">6378</span>
<span style="color: black;"> </span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">##Radians between the xy-plane and the ecliptic plane</span>
<span style="color: black;"> epsilon</span><span style="color: #336366; font-weight: bold;">=</span><span style="color: black; font-weight: bold;">rad(</span><span style="color: purple;">23.45</span><span style="color: black; font-weight: bold;">)</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">##Convert observer's latitude to radians</span>
<span style="color: black;"> L</span><span style="color: #336366; font-weight: bold;">=</span><span style="color: black; font-weight: bold;">rad(</span><span style="color: black;">Lat</span><span style="color: black; font-weight: bold;">)</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## Calculate offset of sunrise based on longitude (min)</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## If Long is negative, then the mod represents degrees West of</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## a standard time meridian, so timing of sunrise and sunset should</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## be made later.</span>
<span style="color: black;"> timezone </span><span style="color: #336366; font-weight: bold;">=</span><span style="color: black;"> </span><span style="color: #803f00;">-</span><span style="color: blue;">4</span><span style="color: #803f00;">*</span><span style="color: black; font-weight: bold;">(abs(</span><span style="color: black;">Long</span><span style="color: black; font-weight: bold;">)</span><span style="color: #803f00;">%%</span><span style="color: blue;">15</span><span style="color: black; font-weight: bold;">)</span><span style="color: #803f00;">*</span><span style="color: black; font-weight: bold;">sign(</span><span style="color: black;">Long</span><span style="color: black; font-weight: bold;">)</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## The earth's mean distance from the sun (km)</span>
<span style="color: black;"> r </span><span style="color: #336366; font-weight: bold;">=</span><span style="color: black;"> </span><span style="color: blue;">149598000</span>
<span style="color: black;"> theta </span><span style="color: #336366; font-weight: bold;">=</span><span style="color: black;"> </span><span style="color: blue;">2</span><span style="color: #803f00;">*</span><span style="color: black;">pi</span><span style="color: #803f00;">/</span><span style="color: purple;">365.25</span><span style="color: #803f00;">*</span><span style="color: black; font-weight: bold;">(</span><span style="color: black;">d</span><span style="color: #803f00;">-</span><span style="color: blue;">80</span><span style="color: black; font-weight: bold;">)</span>
<span style="color: black;"> z.s </span><span style="color: #336366; font-weight: bold;">=</span><span style="color: black;"> r</span><span style="color: #803f00;">*</span><span style="color: black; font-weight: bold;">sin(</span><span style="color: black;">theta</span><span style="color: black; font-weight: bold;">)</span><span style="color: #803f00;">*</span><span style="color: black; font-weight: bold;">sin(</span><span style="color: black;">epsilon</span><span style="color: black; font-weight: bold;">)</span>
<span style="color: black;"> r.p </span><span style="color: #336366; font-weight: bold;">=</span><span style="color: black;"> </span><span style="color: black; font-weight: bold;">sqrt(</span><span style="color: black;">r</span><span style="color: #803f00;">^</span><span style="color: blue;">2</span><span style="color: #803f00;">-</span><span style="color: black;">z.s</span><span style="color: #803f00;">^</span><span style="color: blue;">2</span><span style="color: black; font-weight: bold;">)</span>
<span style="color: black;"> t0 </span><span style="color: #336366; font-weight: bold;">=</span><span style="color: black;"> </span><span style="color: blue;">1440</span><span style="color: #803f00;">/</span><span style="color: black; font-weight: bold;">(</span><span style="color: blue;">2</span><span style="color: #803f00;">*</span><span style="color: black;">pi</span><span style="color: black; font-weight: bold;">)</span><span style="color: #803f00;">*</span><span style="color: black; font-weight: bold;">acos((</span><span style="color: black;">R</span><span style="color: #803f00;">-</span><span style="color: black;">z.s</span><span style="color: #803f00;">*</span><span style="color: black; font-weight: bold;">sin(</span><span style="color: black;">L</span><span style="color: black; font-weight: bold;">))</span><span style="color: #803f00;">/</span><span style="color: black; font-weight: bold;">(</span><span style="color: black;">r.p</span><span style="color: #803f00;">*</span><span style="color: black; font-weight: bold;">cos(</span><span style="color: black;">L</span><span style="color: black; font-weight: bold;">)))</span>
<span style="color: black;"> </span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">##a kludge adjustment for the radius of the sun</span>
<span style="color: black;"> that </span><span style="color: #336366; font-weight: bold;">=</span><span style="color: black;"> t0</span><span style="color: #803f00;">+</span><span style="color: blue;">5</span><span style="color: black;"> </span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## Adjust "noon" for the fact that the earth's orbit is not circular:</span>
<span style="color: black;"> n </span><span style="color: #336366; font-weight: bold;">=</span><span style="color: black;"> </span><span style="color: blue;">720</span><span style="color: #803f00;">-</span><span style="color: blue;">10</span><span style="color: #803f00;">*</span><span style="color: black; font-weight: bold;">sin(</span><span style="color: blue;">4</span><span style="color: #803f00;">*</span><span style="color: black;">pi</span><span style="color: #803f00;">*</span><span style="color: black; font-weight: bold;">(</span><span style="color: black;">d</span><span style="color: #803f00;">-</span><span style="color: blue;">80</span><span style="color: black; font-weight: bold;">)</span><span style="color: #803f00;">/</span><span style="color: purple;">365.25</span><span style="color: black; font-weight: bold;">)</span><span style="color: #803f00;">+</span><span style="color: blue;">8</span><span style="color: #803f00;">*</span><span style="color: black; font-weight: bold;">sin(</span><span style="color: blue;">2</span><span style="color: #803f00;">*</span><span style="color: black;">pi</span><span style="color: #803f00;">*</span><span style="color: black;">d</span><span style="color: #803f00;">/</span><span style="color: purple;">365.25</span><span style="color: black; font-weight: bold;">)</span>
<span style="color: black;"> </span><span style="color: #b00000; font-weight: bold;">## now sunrise and sunset are:</span>
<span style="color: black;"> sunrise </span><span style="color: #336366; font-weight: bold;">=</span><span style="color: black;"> </span><span style="color: black; font-weight: bold;">(</span><span style="color: black;">n</span><span style="color: #803f00;">-</span><span style="color: black;">that</span><span style="color: #803f00;">+</span><span style="color: black;">timezone</span><span style="color: black; font-weight: bold;">)</span><span style="color: #803f00;">/</span><span style="color: blue;">60</span>
<span style="color: black;"> sunset </span><span style="color: #336366; font-weight: bold;">=</span><span style="color: black;"> </span><span style="color: black; font-weight: bold;">(</span><span style="color: black;">n</span><span style="color: #803f00;">+</span><span style="color: black;">that</span><span style="color: #803f00;">+</span><span style="color: black;">timezone</span><span style="color: black; font-weight: bold;">)</span><span style="color: #803f00;">/</span><span style="color: blue;">60</span>
<span style="color: black;"> </span><span style="color: black; font-weight: bold;">return(list(</span><span style="color: #dd0000;">"sunrise"</span><span style="color: black;"> </span><span style="color: #336366; font-weight: bold;">=</span><span style="color: black;"> sunrise,</span><span style="color: #dd0000;">"sunset"</span><span style="color: black;"> </span><span style="color: #336366; font-weight: bold;">=</span><span style="color: black;"> sunset</span><span style="color: black; font-weight: bold;">))</span>
<span style="color: black;">}</span>
</code></pre></pre><span style="font-size: 85%;"><span style="font-family: courier new;"></span></span></div>Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com3tag:blogger.com,1999:blog-1901659770749170850.post-90193858983587155722007-10-14T13:19:00.000-04:002007-10-18T00:35:17.984-04:00Convert polar coordinates to CartesianWhen I want to calculate the coordinates of a location (e.g., a nest or burrow) based on distance and bearing from a grid point, this function helps me avoid writing down SOH-CAH-TOA every time. Just note that the bearing in this case is from the grid point (known location) to the unknown location.<br /><pre><pre><br /><span style="color: rgb(0, 0, 0);">polar2cart</span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="color: rgb(0, 0, 191);">function</span><span style="font-weight: bold; color: rgb(0, 0, 0);">(</span><span style="color: rgb(0, 0, 0);">x,y,dist,bearing,</span><span style="color: rgb(128, 0, 0);">as.deg=</span><span style="color: rgb(0, 128, 0);">FALSE</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">{</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">## Translate Polar coordinates into Cartesian coordinates</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">## based on starting location, distance, and bearing</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">## as.deg indicates if the bearing is in degrees (T) or radians (F)</span><br /><span style="color: rgb(0, 0, 0);"> </span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="color: rgb(0, 0, 191);">if</span><span style="font-weight: bold; color: rgb(0, 0, 0);">(</span><span style="color: rgb(0, 0, 0);">as.deg</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">{</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##if bearing is in degrees, convert to radians</span><br /><span style="color: rgb(0, 0, 0);"> bearing</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="color: rgb(0, 0, 0);">bearing</span><span style="color: rgb(128, 63, 0);">*</span><span style="color: rgb(0, 0, 0);">pi</span><span style="color: rgb(128, 63, 0);">/</span><span style="color: rgb(0, 0, 255);">180</span><br /><span style="color: rgb(0, 0, 0);"> }</span><br /><span style="color: rgb(0, 0, 0);"> </span><br /><span style="color: rgb(0, 0, 0);"> newx</span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="color: rgb(0, 0, 0);">x</span><span style="color: rgb(128, 63, 0);">+</span><span style="color: rgb(0, 0, 0);">dist</span><span style="color: rgb(128, 63, 0);">*</span><span style="font-weight: bold; color: rgb(0, 0, 0);">sin(</span><span style="color: rgb(0, 0, 0);">bearing</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##X</span><br /><span style="color: rgb(0, 0, 0);"> newy</span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="color: rgb(0, 0, 0);">y</span><span style="color: rgb(128, 63, 0);">+</span><span style="color: rgb(0, 0, 0);">dist</span><span style="color: rgb(128, 63, 0);">*</span><span style="font-weight: bold; color: rgb(0, 0, 0);">cos(</span><span style="color: rgb(0, 0, 0);">bearing</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">##Y</span><br /><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">return(list(</span><span style="color: rgb(221, 0, 0);">"x"</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="color: rgb(0, 0, 0);">newx,</span><span style="color: rgb(221, 0, 0);">"y"</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="color: rgb(0, 0, 0);">newy</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><br /><span style="color: rgb(0, 0, 0);">}</span><br /><br /><br /><span style="font-weight: bold; color: rgb(176, 0, 0);">##Example</span><br /><br /><span style="color: rgb(0, 0, 0);">oldloc</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">c(</span><span style="color: rgb(0, 0, 255);">0</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(0, 0, 255);">5</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);"> </span><br /><span style="color: rgb(0, 0, 0);">bearing</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="color: rgb(0, 0, 255);">200</span><span style="color: rgb(0, 0, 0);"> </span><span style="font-style: italic; color: rgb(128, 128, 128);">#degrees</span><br /><span style="color: rgb(0, 0, 0);">dist</span><span style="font-weight: bold; color: rgb(51, 99, 102);">=</span><span style="color: rgb(0, 0, 255);">5</span><br /><br /><span style="color: rgb(0, 0, 0);">newloc</span><span style="font-weight: bold; color: rgb(51, 99, 102);"><-</span><span style="font-weight: bold; color: rgb(0, 0, 0);">polar2cart(</span><span style="color: rgb(0, 0, 0);">oldloc[</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">],oldloc[</span><span style="color: rgb(0, 0, 255);">2</span><span style="color: rgb(0, 0, 0);">],dist,bearing,</span><span style="color: rgb(0, 128, 0);">TRUE</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><span style="font-weight: bold; color: rgb(0, 0, 0);">plot(</span><span style="color: rgb(0, 0, 0);">oldloc[</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(0, 0, 0);">],oldloc[</span><span style="color: rgb(0, 0, 255);">2</span><span style="color: rgb(0, 0, 0);">],</span><span style="color: rgb(128, 0, 0);">xlim=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">c(</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(0, 0, 255);">10</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(0, 0, 255);">10</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(128, 0, 0);">ylim=</span><span style="font-weight: bold; color: rgb(0, 0, 0);">c(</span><span style="color: rgb(128, 63, 0);">-</span><span style="color: rgb(0, 0, 255);">10</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(0, 0, 255);">10</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><br /><span style="font-weight: bold; color: rgb(0, 0, 0);">points(</span><span style="color: rgb(0, 0, 0);">newloc</span><span style="color: rgb(128, 63, 0);">$</span><span style="color: rgb(0, 0, 0);">x,newloc</span><span style="color: rgb(128, 63, 0);">$</span><span style="color: rgb(0, 0, 0);">y,</span><span style="color: rgb(128, 0, 0);">col=</span><span style="color: rgb(221, 0, 0);">"red"</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /></pre></pre><span style="font-size:85%;"><span style="font-family:courier new;"></span></span>Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com3tag:blogger.com,1999:blog-1901659770749170850.post-40539209310110593112007-10-03T16:45:00.000-04:002007-11-29T17:13:13.581-05:00Reorder factor levelsVery often, especially when plotting data, I need to reorder the levels of a factor because the default order is alphabetical. There must be many ways of reordering the levels; however, I always forget which package to look in. A direct way of reordering, using standard syntax is as follows:<br /><div style="overflow:auto"><pre><span style="font-weight: bold; color: rgb(176, 0, 0);">## generate data</span><br /><pre><span style="color: rgb(0, 0, 0);">x </span><span style="font-family:mon;"><span style="font-weight: bold;">=</span></span><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">factor(sample(</span><span style="color: rgb(0, 0, 0);">letters[</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(128, 63, 0);">:</span><span style="color: rgb(0, 0, 255);">5</span><span style="color: rgb(0, 0, 0);">],</span><span style="color: rgb(0, 0, 255);">100</span><span style="color: rgb(0, 0, 0);">, </span><span style="color: rgb(128, 0, 0);">replace=</span><span style="color: rgb(0, 128, 0);">TRUE</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><br /><br /><span style="font-weight: bold; color: rgb(0, 0, 0);">print(levels(</span><span style="color: rgb(0, 0, 0);">x</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">## This will show the levels of x are "Levels: a b c d e"</span><br /><br /><span style="font-weight: bold; color: rgb(176, 0, 0);">## To reorder the levels:<br />## note, if x is not a factor use levels(factor(x))</span><br /><span style="color: rgb(0, 0, 0);">x </span><span style="font-family:mon;"><span style="font-weight: bold;">=</span></span><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(0, 0, 0);">factor(</span><span style="color: rgb(0, 0, 0);">x,</span><span style="font-weight: bold; color: rgb(0, 0, 0);">levels(</span><span style="color: rgb(0, 0, 0);">x</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">[</span><span style="font-weight: bold; color: rgb(0, 0, 0);">c(</span><span style="color: rgb(0, 0, 255);">4</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(0, 0, 255);">5</span><span style="color: rgb(0, 0, 0);">,</span><span style="color: rgb(0, 0, 255);">1</span><span style="color: rgb(128, 63, 0);">:</span><span style="color: rgb(0, 0, 255);">3</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><span style="color: rgb(0, 0, 0);">]</span><span style="font-weight: bold; color: rgb(0, 0, 0);">)</span><br /><br /><span style="font-weight: bold; color: rgb(0, 0, 0);">print(levels(</span><span style="color: rgb(0, 0, 0);">x</span><span style="font-weight: bold; color: rgb(0, 0, 0);">))</span><span style="color: rgb(0, 0, 0);"> </span><span style="font-weight: bold; color: rgb(176, 0, 0);">## Now "Levels: d e a b c"</span><br /></pre></pre></div>In this case, the level order could be set in the first line; however, if there are many levels (and you don't want to type out all of the levels explicitly), the above method is preferred. Again, if there is a better way to do this, please let me know in the comments.Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com1tag:blogger.com,1999:blog-1901659770749170850.post-75210486709885179992007-09-05T10:43:00.000-04:002007-09-05T16:16:08.638-04:00Extract objects from a listWhen using Rmpi to send processes to many nodes, it is convenient to create a list of tasks that are assigned to nodes as they become available. In my case, I was working through a large factorial set of simulations and needed to use a unique set of variable values for each task. Rather than assign these variables individually within the slave function, I decided to extract the variables directly from the list using a variation of the following code:<br /><pre><br /><span style="color: #000000;">mylist</span><span style="font-weight: bold;color: #336366;">=</span><span style="font-weight: bold;color: #000000;">(list(</span><span style="color: #800000;">a=</span><span style="color: #0000ff;">1</span><span style="color: #000000;">,</span><span style="color: #800000;">b=</span><span style="color: #0000ff;">2</span><span style="color: #000000;">,</span><span style="color: #800000;">c=</span><span style="color: #dd0000;">"string1"</span><span style="color: #000000;">,</span><span style="color: #800000;">d=</span><span style="font-weight: bold;color: #000000;">list(</span><span style="color: #dd0000;">"r"</span><span style="font-weight: bold;color: #336366;">=</span><span style="color: #0000ff;">2</span><span style="color: #000000;">,</span><span style="color: #dd0000;">"z"</span><span style="font-weight: bold;color: #336366;">=</span><span style="color: #dd0000;">"string2"</span><span style="font-weight: bold;color: #000000;">)))</span><br /><br /><span style="color: #0000bf;">for</span><span style="font-weight: bold;color: #000000;">(</span><span style="color: #000000;">i </span><span style="color: #0000bf;">in</span><span style="color: #000000;"> </span><span style="color: #0000ff;">1</span><span style="color: #803f00;">:</span><span style="font-weight: bold;color: #000000;">length(</span><span style="color: #000000;">mylist</span><span style="font-weight: bold;color: #000000;">))</span><span style="color: #000000;">{</span><br /><span style="color: #000000;"> </span><span style="font-weight: bold;color: #b00000;">##first extract the object value</span><br /><span style="color: #000000;"> tempobj</span><span style="font-weight: bold;color: #336366;">=</span><span style="color: #000000;">mylist[[i]]</span><br /><span style="color: #000000;"> </span><span style="font-weight: bold;color: #b00000;">##now create a new variable with the original name of the list item</span><br /><span style="color: #000000;"> </span><span style="font-weight: bold;color: #000000;">eval(parse(</span><span style="color: #800000;">text=</span><span style="font-weight: bold;color: #000000;">paste(names(</span><span style="color: #000000;">mylist</span><span style="font-weight: bold;color: #000000;">)</span><span style="color: #000000;">[[i]],</span><span style="color: #dd0000;">"= tempobj"</span><span style="font-weight: bold;color: #000000;">)))</span><br /><span style="color: #000000;">}</span><br /><br /><span style="color: #000000;">> </span><span style="font-weight: bold;color: #000000;">print(</span><span style="color: #000000;">a</span><span style="font-weight: bold;color: #000000;">)</span><br /><span style="color: #000000;">[</span><span style="color: #0000ff;">1</span><span style="color: #000000;">] </span><span style="color: #0000ff;">1</span><br /><span style="color: #000000;">> </span><span style="font-weight: bold;color: #000000;">print(</span><span style="color: #000000;">b</span><span style="font-weight: bold;color: #000000;">)</span><br /><span style="color: #000000;">[</span><span style="color: #0000ff;">1</span><span style="color: #000000;">] </span><span style="color: #0000ff;">2</span><br /><span style="color: #000000;">> </span><span style="font-weight: bold;color: #000000;">print(</span><span style="color: #000000;">c</span><span style="font-weight: bold;color: #000000;">)</span><br /><span style="color: #000000;">[</span><span style="color: #0000ff;">1</span><span style="color: #000000;">] </span><span style="color: #dd0000;">"string1"</span><br /><span style="color: #000000;">> </span><span style="font-weight: bold;color: #000000;">print(</span><span style="color: #000000;">d</span><span style="font-weight: bold;color: #000000;">)</span><br /><span style="color: #803f00;">$</span><span style="color: #000000;">r</span><br /><span style="color: #000000;">[</span><span style="color: #0000ff;">1</span><span style="color: #000000;">] </span><span style="color: #0000ff;">2</span><br /><br /><span style="color: #803f00;">$</span><span style="color: #000000;">z</span><br /><span style="color: #000000;">[</span><span style="color: #0000ff;">1</span><span style="color: #000000;">] </span><span style="color: #dd0000;">"string2"</span><br /></pre><br />If there is a more elegant way to do this, please let me know in the comments.Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com2tag:blogger.com,1999:blog-1901659770749170850.post-64359344969118001002007-09-05T10:26:00.001-04:002009-01-20T11:03:39.333-05:00Parallel computing in RLately I have been running R in parallel on a Beowulf cluster using the <a href="http://www.stats.uwo.ca/faculty/yu/Rmpi">Rmpi</a> package. Instrumental in my ability to do this were two excellent tutorials:<br /><br /><a href="http://www.biostat.ucsf.edu/biostat/sen/cluster/RMPI-cluster2.html">The Division of Biostatistics at UCSF</a> provides some of the background required to get things running on the Linux side of things. (Note: sometimes the url for this site comes up with an extra "biostat" in the address -- just delete that and the page should load.<br /><br /><a href="http://math.acadiau.ca/ACMMaC/Rmpi/index.html">The Acadia Centre for Mathematical Modelling and Computation</a> has a variety of sample scripts to point you in the right direction. Once you read through the code, it becomes very easy to set up a template you can re-use in the future.Foresterhttp://www.blogger.com/profile/00291122844817714304noreply@blogger.com0