[R] CIE Diagram with wavelengths tick marks

2023-07-20 Thread Kenneth Knoblauch
Hi,

I'm assuming that by CIE Diagram, you mean the 1931 2 deg standard.  Tables of 
these (and other CIE standards) can be found at the site http://www.cvrl.org 
maintained by Andrew Stockman.  You can find and download these classic CIE 
data and more modern versions from there.  

Since you have plotted the spectrum locus of the chromaticity diagram (I'm 
supposing), you ought to be able to just pick off the wavelengths that you want 
and add a small point or a small line as a tick mark.  If you want it oriented, 
you'll have to do a little calculation to get the direction perpendicular to 
the spectrum locus at each wavelength, but that's not too complicated.

I recently needed a figure for a talk and downloaded the 2 deg 1931 data from 
the above website and used points to indicate specific wavelengths that 
interested me.  The code is below:

d <- read.csv("ciexyz31_1.csv", header = FALSE, stringsAsFactors = TRUE)
names(d) <- c("Wavelength", "X", "Y", "Z")
d$x <- with(d, X/(X + Y + Z))
d$y <-  with(d, Y/(X + Y + Z))

par(pty = "s")
plot(y ~ x, d, type = "l", axes = 'FALSE',
xlim = c(0, 1), ylim = c(0, 1))
lines(y[c(1, nrow(d))]~ x[c(1, nrow(d))], d)
axis(1, seq(0, 1, 0.2))
axis(2, seq(0, 1, 0.2))
lines(c(0, 1), c(1, 0), col = "blue", lwd = 2, lty = 2)

points(d[d$Wavelength == 555, ]$x, d[d$Wavelength == 555, ]$y, pch = 16, col = 
"green")
points(d[d$Wavelength == 589, ]$x, d[d$Wavelength == 589, ]$y, pch = 16, col = 
"yellow")
points(d[d$Wavelength == 670, ]$x, d[d$Wavelength == 670, ]$y, pch = 16, col = 
"red")

Good luck.

best,

Ken


# Hy,

# for plotting CIE Diagrams i found Package pavo with function cieplot. 
# That works fine.
# Now i want to have wavelength tick marks as well around the plotting area.
# Is there a way to do so, also other ways/Pakages to plot are welcome.

# My sample Code:
   # library(pavo)

   # coldat2 <- as.data.frame(matrix(rep(1/3, 3), nrow = 1, ncol = 3))

   # # Make sure this dataset works with the cieplot() function
   # attr(coldat2, "clrsp") <- "CIEXYZ"
   # colnames(coldat2) <- c("x", "y", "z")

   # cieplot(coldat2, col="white", main="CIE Test Plot")



# Best regards
# Tilmann



 
  
 
  ___
  Kenneth Knoblauch
 Inserm U1208
 Stem-cell and Brain Research Institute
 18 avenue du Doyen Lépine
 69500 Bron
 France
 tel: +33 (0)4 72 91 34 77
 fax: +33 (0)4 72 91 34 61
 portable: +33 (0)6 84 10 64 10
 https://sbri.fr/public-profile/63/single-member/

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] lines through points in lattice legend

2023-01-28 Thread Kenneth Knoblauch


Thanks.  That works well and it's simple.  

Ken

  ___
  Kenneth Knoblauch
 Inserm U1208
 Stem-cell and Brain Research Institute
 18 avenue du Doyen Lépine
 69500 Bron
 France
 tel: +33 (0)4 72 91 34 77
 fax: +33 (0)4 72 91 34 61
 portable: +33 (0)6 84 10 64 10
 https://sbri.fr/public-profile/63/single-member/





De : Deepayan Sarkar 
Envoyé : samedi 28 janvier 2023 16:46
À : Kenneth Knoblauch
Cc : r-help@r-project.org
Objet : Re: [R] lines through points in lattice legend
    
On Sat, Jan 28, 2023 at 2:49 PM Kenneth Knoblauch
 wrote:
>
> Hi,
>
> I'm struggling to find if there is a simple way to make the lines and points 
> overlap in a legend for a lattice plot using auto.key.  Here is a toy example 
> of what doesn't work (for me) as the lines and points are adjacent rather 
> than overlapping:
>
> library(lattice)
> d <- data.frame(x = 1:2, y = 1:4, f = factor(letters[1:2]))
>
> xyplot(y ~ x, d, groups = f, type = "b",
> pch = 21:22, fill = c("white", "black"), col = "black",
> par.settings = list(superpose.symbol =
> list(pch = 21:22, fill = c("white", "black"), col = "black"),
> superpose.line = list(col = "black")),
> auto.key = list(space = "right", points = TRUE, lines = TRUE))

Just adding a type = "b" should work (for the lines), so something like:

   auto.key = list(space = "right", type = "b",
   points = FALSE, lines = TRUE)

BTW, once you specify par.settings, you don't need to specify the
parameters again, so you can drop the second line.

Best,
-Deepayan

>
>
> I've seen a few examples on stack.overflow but I haven't been able to make 
> them work, or they seem more complicated then I would think it should be (but 
> I don't exclude that I'm fooling myself there).
>
> R version 4.2.2 Patched (2023-01-04 r83564)
> Platform: x86_64-apple-darwin17.0 (64-bit)
> Running under: macOS Catalina 10.15.7
>
> Matrix products: default
> BLAS:   
> /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
> LAPACK: 
> /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
>
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>
> attached base packages:
> [1] stats graphics  grDevices utils datasets  methods   base
>
> other attached packages:
> [1] lattice_0.20-45
>
> loaded via a namespace (and not attached):
> [1] compiler_4.2.2 grid_4.2.2
>
>  Thanks for any enlightenment, in advance.
>
> Ken
>
>   ___
>   Kenneth Knoblauch
>  Inserm U1208
>  Stem-cell and Brain Research Institute
>  18 avenue du Doyen Lépine
>  69500 Bron
>  France
>  tel: +33 (0)4 72 91 34 77
>  fax: +33 (0)4 72 91 34 61
>  portable: +33 (0)6 84 10 64 10
>  https://sbri.fr/public-profile/63/single-member/
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide  http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] lines through points in lattice legend

2023-01-28 Thread Kenneth Knoblauch
Hi,

I'm struggling to find if there is a simple way to make the lines and points 
overlap in a legend for a lattice plot using auto.key.  Here is a toy example 
of what doesn't work (for me) as the lines and points are adjacent rather than 
overlapping:

library(lattice)
d <- data.frame(x = 1:2, y = 1:4, f = factor(letters[1:2]))

xyplot(y ~ x, d, groups = f, type = "b", 
pch = 21:22, fill = c("white", "black"), col = "black",
par.settings = list(superpose.symbol = 
list(pch = 21:22, fill = c("white", "black"), col = "black"),
superpose.line = list(col = "black")),
auto.key = list(space = "right", points = TRUE, lines = TRUE))

I've seen a few examples on stack.overflow but I haven't been able to make them 
work, or they seem more complicated then I would think it should be (but I 
don't exclude that I'm fooling myself there).

R version 4.2.2 Patched (2023-01-04 r83564)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   
/Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
LAPACK: 
/Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base 

other attached packages:
[1] lattice_0.20-45

loaded via a namespace (and not attached):
[1] compiler_4.2.2 grid_4.2.2
 
 Thanks for any enlightenment, in advance.

Ken 
 
  ___
  Kenneth Knoblauch
 Inserm U1208
 Stem-cell and Brain Research Institute
 18 avenue du Doyen Lépine
 69500 Bron
 France
 tel: +33 (0)4 72 91 34 77
 fax: +33 (0)4 72 91 34 61
 portable: +33 (0)6 84 10 64 10
 https://sbri.fr/public-profile/63/single-member/

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] Need Help To Solve An Equation In R

2017-05-28 Thread Kenneth Knoblauch
This looks an awful lot like you are trying to solve for d' in an 
m-alternative forced choice experiment for an unbiased observer.  Try 
the function dprime.mAFC from the psyphy package.  For comparison with 
your example, it's code is:


 psyphy:::dprime.mAFC
function (Pc, m)
{
m <- as.integer(m)
if (m < 2)
stop("m must be an integer greater than 1")
if (!is.integer(m))
stop("m must be an integer")
if (Pc <= 0 || Pc >= 1)
stop("Pc must be in (0,1)")
est.dp <- function(dp) {
pr <- function(x) dnorm(x - dp) * pnorm(x)^(m - 1)
Pc - integrate(pr, lower = -Inf, upper = Inf)$value
}
dp.res <- uniroot(est.dp, interval = c(-10, 10))
dp.res$root
}


Hope that helps,

best,

Ken


On Sat, May 27, 2017 at 9:16 AM, Neetu Shah  
wrote:

Dear Sir/Ma'am,

I am trying to make a function to solve an equation that is given 
below:

findH <- function(p_star, k){
fun <- function(y){
(pnorm(y+h))^(k-1)*dnorm(y)
}
lhs <- integrate(fun, -Inf, Inf)$value
return(uniroot(lhs, lower = -10, upper = 10,
tol = p_star)$root)
}
In lhs, I integrated a function with respect to y and there would be an
unknown h that I need to find.
I need to equate this lhs to p_star value in order to find h. Which
function should I use to solve this equation?
Please guide me regarding this.
Thank you.

--
With Regards,
Neetu Shah,
B.Tech.(CS),
3rd Year,
IIIT Vadodara.

--
Kenneth Knoblauch
Inserm U1208
Stem-cell and Brain Research Institute
18 avenue du Doyen Lépine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.sbri.fr/members/kenneth-knoblauch.html

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] meaning of lm( y~., data=mydat ), is it a language feature, is it documented, is it supported?

2016-05-23 Thread Kenneth Knoblauch
John Sorkin  grecc.umaryland.edu> writes:
> The syntax
> mydat <- data.frame( y,x )
> fit1 <- lm( y~., data=mydat )
> appears to perform a multivariable regression of y on 
every non-y variable in the data frame mydat. I can not
> find this syntax (y~.) in R documentation. Is y~. 
a supported feature of the R language? Where can I find it
> documented? I would hate to write code that
 is dependent on a non-supported, non-documented 
language feature.
> Thank you,
> John
> John David Sorkin M.D., Ph.D.
> Professor of Medicine
> Chief, Biostatistics and Informatics
> University of Maryland School of 
Medicine Division of Gerontology and Geriatric Medicine
> Baltimore VA Medical Center
> 10 North Greene Street
> GRECC (BT/18/GR)
> Baltimore, MD 21201-1524
> (Phone) 410-605-7119
> (Fax) 410-605-7913
> 

How about section 11.5 of An Introduction to R?

-- 
Kenneth Knoblauch
Inserm U1208
Stem-cell and Brain Research Institute
18 avenue du Doyen Lépine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.sbri.fr/members/kenneth-knoblauch.html
__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] speeding up sum of squared differences calculation

2013-10-22 Thread Kenneth Knoblauch
Actually, you don't need to use c() in the dist example, either, which 
should shave off a few microseconds for a function call.  It was late 
last night ...


Ken

On 22-10-2013 15:08, S Ellison wrote:

Conclusion: hard to beat outer() for this purpose in R

... unless you use Ken Knoblauch's suggestion of dist() or Rccp.

Nice indeed.


S Ellison


***
This email and any attachments are confidential. Any use, copying or
disclosure other than by the intended recipient is unauthorised. If
you have received this message in error, please notify the sender
immediately via +44(0)20 8943 7000 or notify postmas...@lgcgroup.com
and delete this message and any copies from your computer and network.
LGC Limited. Registered in England 2991879.
Registered office: Queens Road, Teddington, Middlesex, TW11 0LY, UK


--
Kenneth Knoblauch
Inserm U846
Stem-cell and Brain Research Institute
Department of Integrative Neurosciences
18 avenue du Doyen Lépine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.sbri.fr/members/kenneth-knoblauch.html

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] difference between assignment syntax - vs =

2009-02-23 Thread Kenneth Knoblauch

It's easier to read.  Better machine-human interaction.

ergonomic:  (esp. of workplace design) intended to provide optimum  
comfort and to avoid stress or injury.


Quoting Wacek Kusnierczyk waclaw.marcin.kusnierc...@idi.ntnu.no:


Ken Knoblauch wrote:

Wacek Kusnierczyk Waclaw.Marcin.Kusnierczyk at idi.ntnu.no writes:



Thomas Lumley wrote:


Although it's probably true that most long-time R users use -, this
is at least in part because a long-time R user would initially have
had to use -, since = wasn't available in the distant past.

I would say that it's entirely a matter of taste -- the things that
otherwise could have been traps are mostly syntax errors.  The only
proviso is that if you post code using = it is (even more) important
to leave spaces around the = than it would be for -.


the reason being ...?

vQ




ergonomy!




ergonomy?

vQ





--
Ken Knoblauch
Inserm U846
Stem-cell and Brain Research Institute
Department of Integrative Neurosciences
18 avenue du Doyen Lépine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.sbri.fr/members/kenneth-knoblauch.html

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] length 1 offset in glm

2009-02-23 Thread Kenneth Knoblauch

Hi,

I'm trying to use an offset with glm.  According to the glm man
page

offset	... This should be NULL or a numeric vector of length either  
one or equal to the number of cases. ...


but with the following example, I get an error if the offset is of length 1

c1 - structure(list(Contr = c(0.028, 0.043, 0.064, 0.097, 0.146, 0.219
+ ), Correct = c(34L, 57L, 94L, 152L, 160L, 160L), Incorrect = c(126L,
+ 103L, 66L, 8L, 0L, 0L)), .Names = c(Contr, Correct, Incorrect
+ ), row.names = c(NA, 6L), class = data.frame)

glm(cbind(Correct, Incorrect) ~ Contr - 1, binomial,
data = c1, offset = qlogis(0.25))

Error in model.frame.default(formula = cbind(Correct, Incorrect) ~ Contr -  :
  variable lengths differ (found for '(offset)')

but not if it has the length of the number of rows of the data frame

glm(cbind(Correct, Incorrect) ~ Contr - 1, binomial,
data = c1, offset = rep(qlogis(0.25), nrow(c1)))

...
Coefficients:
Contr
28.34

Degrees of Freedom: 6 Total (i.e. Null);  5 Residual
Null Deviance:  1342
Residual Deviance: 92.35AIC: 114.3

(I'm not pushing this as a good model...)

I traced the error to this line in glm

mf - eval(mf, parent.frame())

and the error can be fixed if the following code from later in glm is put
before the above line with the additional line that I added to
(re)set the offset of mf

if (!is.null(offset)) {
if (length(offset) == 1)
offset - rep(offset, NROW(Y))
else if (length(offset) != NROW(Y))
stop(gettextf(number of offsets is %d should equal %d  
(number of observations),

length(offset), NROW(Y)), domain = NA)
}

 mf$offset - offset

I don't know if this breaks something important or is a terrible
hack.  Or, am I misreading the way glm should work...

sessionInfo()
R version 2.8.1 Patched (2009-02-22 r47985)
i386-apple-darwin8.11.1

locale:
en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

loaded via a namespace (and not attached):
[1] tools_2.8.1

Thanks for any enlightenment.

Ken


--
Ken Knoblauch
Inserm U846
Stem-cell and Brain Research Institute
Department of Integrative Neurosciences
18 avenue du Doyen Lépine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.sbri.fr/members/kenneth-knoblauch.html

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


[R] puzzle about contrasts

2008-09-09 Thread Kenneth Knoblauch

Hi,

I'm trying to redefine the contrasts for a linear model.
With a 2 level factor, x, with levels A and B, a two level
factor outputs A and B - A from an lm fit, say
lm(y ~ x). I would like to set the contrasts so that
the coefficients output are -0.5 (A + B) and B - A,
but I can't get the sign correct for the first coefficient
(Intercept).

Here is a toy example,

set.seed(12161952)
y - rnorm(10)
x - factor(rep(letters[1:2], each = 5))
##  so  A and B =
tapply(y, x, mean)

a  b
-0.719  0.8323837

## and with treatment contrasts
coef(lm(y ~ x))  ## A and B - A

(Intercept)  xb
 -0.719   1.5522724

Then, I try to redefine the contrasts

### would like contrasts: -0.5 (A + B) and B - A
D1 - matrix( c(-0.5, -0.5,
-1, 1),
2, 2, byrow = TRUE)
C1 - solve(D1)
Cnt - C1[, -1]
contrasts(x) - Cnt
coef(lm(y ~ x))

(Intercept)  x1
 0.05624745  1.55227241

but note that the desired value is
-0.5 * sum(tapply(y, x, mean))

[1] -0.05624745

I note that the first column of C1 is -1's not +1's
and that working by hand, if I tamper with the model matrix

mm - model.matrix(y ~ x)
mm[, 1] - -1

mm
   (Intercept)   x1
1   -1 -0.5
2   -1 -0.5
3   -1 -0.5
4   -1 -0.5
5   -1 -0.5
6   -1  0.5
7   -1  0.5
8   -1  0.5
9   -1  0.5
10  -1  0.5
attr(,assign)
[1] 0 1
attr(,contrasts)
attr(,contrasts)$x
  [,1]
a -0.5
b  0.5

solve(t(mm) %*% mm) %*% t(mm) %*% y  ##Yes, I know. Use QR
   [,1]
(Intercept) -0.05624745
x1   1.55227241

gives the correct sign.

So, I guess my question reduces to how one would set the
contrasts for the model.matrix to be correct
for this to work out correctly?

Thank you.

Ken


--
Ken Knoblauch
Inserm U846
Institut Cellule Souche et Cerveau
Département Neurosciences Intégratives
18 avenue du Doyen Lépine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.sbri.fr

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.