Custom function applied to elements in matrix

jnr20

New Member
#1
So i have the function: p^x/(1-p)^y.
I want to produce a 10x10 matrix, where each element within the matrix is evaluated by the above function. where x and y=1,2,..,10, and p=0.3 for example.
So the matrix at (1,1) would be (0.3^1 / (1-0.3)^1)=0.4285714
at (1,2) would be (0.3^1 / (1-0.3)^2)=0.6122449
and so on.

Code:
x=c(1,2,3,4,5,6,7,8,9,10)
y=c(1,2,3,4,5,6,7,8,9,10)
func<-function(p,x,y){p^x / (1-p)^y}
matrix(data=func(0.3,x,y),ncol=10,nrow=10)
the code above produces the following (only first 2 cols and rows shown):
Code:
          [,1]           [,2]       
 [1,] 0.4285714286 0.4285714286
 [2,] 0.1836734694 0.1836734694
 [3,] 0.0787172012 0.0787172012
 [4,] 0.0337359434 0.0337359434 
 [5,] 0.0144582614 0.0144582614
this isn't what i'm after.
i tried to use the function outer, but got an error.
Code:
> outer(x,y,FUN='func(0.3,x,y)')
Error in get(as.character(FUN), mode = "function", envir = envir) : 
  object 'func(0.3,x,y)' of mode 'function' was not found
what am i doing wrong? and how do i get the desired results?

thanks :)
 

Jake

Cookie Scientist
#2
First create the matrix of x and y input values, and create a scalar for p, then just call your function directly on x, y, and p. This should return a matrix of results.
 

bryangoodrich

Probably A Mammal
#3
It seems to me more efficient to think of the data in "long format" where each row is the given pair you're interested in. Then you can apply your function to each pair, appropriately.

Code:
x <- 1:10
y <- 1:10
z <- expand.grid(x, y)  # creates all the pairs
apply(z, 1, function(pair, p) {p^pair[1] / (1-p)^pair[2]}, p = 0.3)  # Since z has 2 columns, each row ('pair') is a 2-vector (x,y)
This is certainly not the "matrix" approach you were interested in, but it seems the most intuitive R way to go about it. Columns should reflect unique information about each row. A matrix as you want to use it reflects individual pieces of information determined by the row and column (wide format). With the 'z' object I created, you'll see now that each row is that same piece of information and each column represents that unique piece of information (row and column) that the matrix dimensions provided. Having data in this long format, then, offers a benefit for information processing. In this case, you can simply supply your function across the rows to get what you want. If you change your function to something like

Code:
MyFun <- function(pair, p)
{
    x <- pair[1]
    y <- pair[2]
    p^x / (1-p)^y
}
You can simplify my above code to

Code:
apply(z, 1, MyFunc, p=0.3)
Because now the function handles the 2-tuple (vector) appropriately for the process you're interested in.

Update: To get the result as a matrix as in Jake's method, you simply need to convert the result

Code:
round(matrix(apply(z,1,MyFunc), ncol=10), 4)  # Rounding optional
#         [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]   [,10]
#  [1,] 0.4286 0.6122 0.8746 1.2495 1.7850 2.5500 3.6428 5.2040 7.4343 10.6204
#  [2,] 0.1286 0.1837 0.2624 0.3748 0.5355 0.7650 1.0928 1.5612 2.2303  3.1861
#  [3,] 0.0386 0.0551 0.0787 0.1125 0.1606 0.2295 0.3279 0.4684 0.6691  0.9558
#  [4,] 0.0116 0.0165 0.0236 0.0337 0.0482 0.0688 0.0984 0.1405 0.2007  0.2868
#  [5,] 0.0035 0.0050 0.0071 0.0101 0.0145 0.0207 0.0295 0.0422 0.0602  0.0860
#  [6,] 0.0010 0.0015 0.0021 0.0030 0.0043 0.0062 0.0089 0.0126 0.0181  0.0258
#  [7,] 0.0003 0.0004 0.0006 0.0009 0.0013 0.0019 0.0027 0.0038 0.0054  0.0077
#  [8,] 0.0001 0.0001 0.0002 0.0003 0.0004 0.0006 0.0008 0.0011 0.0016  0.0023
#  [9,] 0.0000 0.0000 0.0001 0.0001 0.0001 0.0002 0.0002 0.0003 0.0005  0.0007
# [10,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0001 0.0001  0.0002
 

bryangoodrich

Probably A Mammal
#4
To use your function, you'd also want to define it

Code:
f <- function(x,y, p=0.3)
and use outer like my apply above

Code:
outer(x, y, f, p=0.3)  # or leave p out in this case if 0.3 is the default you like
 

Jake

Cookie Scientist
#5
Here's the vectorized method I hinted at:
Code:
x <- matrix(1:10, nrow=10, ncol=10)
y <- matrix(1:10, nrow=10, ncol=10, byrow=T)
p <- .3
round(p^x/(1-p)^y, digits=4)
#         [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]   [,10]
#  [1,] 0.4286 0.6122 0.8746 1.2495 1.7850 2.5500 3.6428 5.2040 7.4343 10.6204
#  [2,] 0.1286 0.1837 0.2624 0.3748 0.5355 0.7650 1.0928 1.5612 2.2303  3.1861
#  [3,] 0.0386 0.0551 0.0787 0.1125 0.1606 0.2295 0.3279 0.4684 0.6691  0.9558
#  [4,] 0.0116 0.0165 0.0236 0.0337 0.0482 0.0688 0.0984 0.1405 0.2007  0.2868
#  [5,] 0.0035 0.0050 0.0071 0.0101 0.0145 0.0207 0.0295 0.0422 0.0602  0.0860
#  [6,] 0.0010 0.0015 0.0021 0.0030 0.0043 0.0062 0.0089 0.0126 0.0181  0.0258
#  [7,] 0.0003 0.0004 0.0006 0.0009 0.0013 0.0019 0.0027 0.0038 0.0054  0.0077
#  [8,] 0.0001 0.0001 0.0002 0.0003 0.0004 0.0006 0.0008 0.0011 0.0016  0.0023
#  [9,] 0.0000 0.0000 0.0001 0.0001 0.0001 0.0002 0.0002 0.0003 0.0005  0.0007
# [10,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0001 0.0001  0.0002