Solving a system of integral equations in R using multiroot
I have a system of integrals that have free variables S1 and S2. I would like to know what values of S1, S2 make the resulting (integrals1) equal to 0. (I will eventually change the parameter values so that they are not the same equations.)
Multiroot uses the library: library(rootSolve)
.
a11 < 1 #alpha_{11}
a12 < 1 #alpha_{12}
a21 < 1 #alpha_{21}
a22 < 1 #alpha_{22}
b1 < 2 #beta1
b2 < 2 #beta2
d1 < 1 #delta1
d2 < 1 #delta2
g < 0.5 #gamma
integrand1 < function(x,S) {b1*g/d1*exp(g*x)*(1exp(d1* x))*exp(a11*b1*S[1]/d1*(1exp(d1*x))a12*b2*S[2]/d2*(1exp(d2*x)))}
integrand2 < function(x,S) {b2*g/d2*exp(g*x)*(1exp(d2* x))*exp(a22*b2*S[2]/d2*(1exp(d2*x))a21*b1*S[1]/d1*(1exp(d1*x)))}
#defining equation we would like to solve
intfun1< function(S) {integrate(function(x)integrand1(x, S),lower=0,upper=Inf)[[1]]1}
intfun2< function(S) {integrate(function(x)integrand2(x, S),lower=0,upper=Inf)[[1]]1}
#putting both equations into one term
model < function(S) c(F1 = intfun1,F2 = intfun2)
#Solving for roots
(ss <multiroot(f=model, start=c(0,0)))
Gives the error:
Error in stode(y, times, func, parms = parms, ...) : REAL() can only be applied to a 'numeric', not a 'list'
For reference this simpler system works fine:
#Defining the functions
model < function(x) c(F1 = x[1]+ 4*x[2] 8,F2 = x[1]4*x[2])
#Solving for the roots
(ss < multiroot(f = model, start = c(0,0)))
If I am mistaken in using multiroot for solving the roots of a system of equations and you think there is a better function I am open to suggestions. Thank you!
See also questions close to this topic

R  get values of matrix using filter in other matrix
I have a matrix of the following type:
c1 c2 c3 c4 1234 123 124 134 234 23 14 124 1234 123 124 134 234 23 14 124 I have another matrix that tells me what are the cells that are ok:
c1 c2 c3 c4 1 1 1 1 1 0 1 1 0 0 1 1 0 1 0 1 I want to get values of the 2 first columnes on table 1, that have a 1 on table 2:
c1 c2 1234 123 234 14 124 134 23 124 I can solve it with a for lop and head, but is there a better way?

Copy values from one row to another based on condition
I have this dataset:
# Data movmnt_id < c("101", "601", "105", "321") plant < c("FF", "FF", "DO", "BO") loc < c("MM", "MM", "KB", "RD") vendor < c(123, NA,NA,NA) customer < c(456, NA,NA,NA) check < c(NA, NA, "defined", "defined") df < data.frame(movmnt_id, plant, loc, vendor,customer,check) movmnt_id plant loc vendor customer check 1 101 FF MM 123 456 <NA> 2 601 FF MM NA NA <NA> 3 105 DO KB NA NA defined 4 321 BO RD NA NA defined
I need to get this output (in the second row
vendor
&customer
are copied from the first row):movmnt_id plant loc vendor customer check 1 101 FF MM 123 456 <NA> 2 601 FF MM 123 456 <NA> 3 105 DO KB NA NA defined 4 321 BO RD NA NA defined
The condition is following:
If in current row `movmnt_id `== 601 > take row *WHERE* `plant` & `loc` are the same as in the current row *AND* `movmnt_id == 101` *AND* is.na(check) > copy from found row `vendor` & `customer` to the current row
I could think about some for loop, but for my dataset it will be too heavy.
I am wondering if there is more elegant solution with less computational cost.I was trying to adapt solutions from these cases, but did not succeed:

How to Edit Code to Query Google Street View API for Entire City?
Using the following code, you can query the Google Street View API which returns a single image from the coordinates specified (
51.509865, 0.118092
):# install.packages('googleway') library(googleway) api_key < 'ENTER API KEY' google_streetview( location = c(51.509865, 0.118092), # lat/lon coordinates size = c(600, 400), # w x h (image size) heading = 0, # north facing camera pitch = 0, # camera angle output = 'plot', response_check = TRUE, key = api_key )
Documentation: https://cran.rproject.org/web/packages/googleway/googleway.pdf
However, this returns a single image from the coordinates specified. How could I amend the code so that it returns 10 images per square kilometre of an entire city? For instance, London (UK) is 1,572 km²  I want to return 10 images per square kilometre so that I have 15,720 images covering the whole of London.
I suspect a loop might be required, but I'm not sure how it would work in practise because you'd need to adjust the coordinates with each loop. Any ideas on how to adjust the coordinates so that it stays within the boundaries of London?

number of items to replace is not a multiple of , MonteCarlo Stochastic Process R error
## set.seed(123) SimpleEulerApproximation = function(T,x,a,b,delta){ numberofSteps = T/delta; TimeSteps = rep(numberofSteps,1); Y = rep(numberofSteps,1) Y[1] = x; for (i in 1:numberofSteps){ TimeSteps[i] = 0 + i*delta; } for (j in 2:numberofSteps){ Y[j] = Y[j1] + a*Y[j1]*delta + b*Y[j1]*rnorm(1,0,sqrt(delta)); } ##plot(TimeSteps,Y, type = "l") } SimpleEulerApproximation(1,20,0.01,0.25,0.001) set.seed(123) MultipleEulerApproximation = function(T,x,a,b,delta,numberofTrajectories){ numberofSteps = round(T/delta); TimeSteps = rep(numberofSteps,1); Y = rep(numberofSteps,rep(numberofTrajectories)) Y = data.matrix(Y) for (i in 1:numberofTrajectories){ Y[,i] = SimpleEulerApproximation(T,x,a,b,delta); } for (i in 1:numberofSteps){ TimeSteps[i] = 0 + i*delta; } AverageTrajectory = rep(numberofSteps,1) for (i in 1:numberofSteps){ AverageTrajectory[i] = mean(Y[i,]) } ##plot(TimeSteps,AverageTrajectory) } MultipleEulerApproximation(1,52,0.12,0.30,0.0001,10000) MonteCarloSimulation = function(T,x,r,sigma,K,delta,numberofTrajectories){ Y = MultipleEulerApproximation(T,x,r,sigma,delta,numberofTrajectories); lastStep = round(T/delta); max(Y[lastStep,]K,0); size(Y) price = 1/numberofTrajectories * sum(max(Y[lastStep,]K,0))*exp(r*T) } MonteCarloSimulation(0.25,52,0.12,0.3,50,0.0001,10000)
When I run the code for multipleEulerApproximation, I get replacement has length 0 error. Can someone help me with this? Much Appreciated. The first one is simple Euler Approximation for stochastic differential equation dXt = −0.1Xtdt + 0.25XtdBt, X0 = 20 over the time interval [0, 1] with time step size ∆ = 0.001. The second chunk of code is for multipleeulerapproximation that is where the error. The thirdchunk is for calculating European call option price using projections.

Composite Trapezoidal Rule and the number of nodes needed for the Standard Normal Distribution
I am using different numerical methods to approximate the normal distribution. My first method is Trapezoid Rule. Why is my number of nodes required behaving so oddly? I used a table to find the lowest value for n so that the error would be less than .0001. I would have expectedn a pattern. Is there one and I am missing it? n for a standard deviation of 1 was 6. n for a standard deviation of 2 was 54. n for a standard deviation of 3 was 29.
My code follows:
s=1 print("For Standard Deviation, s, of 1, we have a=1 and b=1") def CompositeTrapezoidal(f, a, b, n): #compute the size of each subinterval h = (ba)/n #here we add all the terms sum = 1/2*f(a) #add the values for all intermediate nodes for i in [1,2,..,n1] : sum += f(a+i*h) #add the last node sum +=1/2*f(b) #multiply by h to complete the computation sum = sum*h #return the sum return sum #define the function we will be using f(x) =(e^((x^2)/2))/(sqrt(2*pi)) #print the top part of the table print ("n", "\t","trapezoidal" , "\t\t", "error") print ("") #call the function several times for power in [11,12,..,45]: n = power trapezoidal = CompositeTrapezoidal(f, 1, 1, n).n() error = (trapezoidal  integral(f(x), x,1,1)).n() if abs(error) >= 0.0001: print (n, "\t", trapezoidal, "\t", error) if abs(error) <= 0.0001: print(n, "nodes") print("by the table above we need n=41 for s=1 to get error less than 1/100.") trapezoidal = CompositeTrapezoidal(f, 1, 1, 41).n() error = (trapezoidal  integral(f(x), x,1,1)).n() print("integral(f, x, 1, 1) ≈", trapezoidal, "by Composite Trapezoid Rule") print() #s=2 print("For Standard Deviation, s, of 2, we have a=2 and b=2") def CompositeTrapezoidal(f, a, b, n): #compute the size of each subinterval h = (ba)/n #here we add all the terms sum = 1/2*f(a) #add the values for all intermediate nodes for i in [1,2,..,n1] : sum += f(a+i*h) #add the last node sum +=1/2*f(b) #multiply by h to complete the computation sum = sum*h #return the sum return sum #define the function we will be using f(x) =(e^((x^2)/2))/(sqrt(2*pi)) #print the top part of the table print ("n", "\t","trapezoidal" , "\t\t", "error") print ("") #call the function several times for power in [21,22,..,55]: n = power trapezoidal = CompositeTrapezoidal(f, 2, 2, n).n() error = (trapezoidal  integral(f(x), x,2,2)).n() if abs(error) >= 0.0001: print (n, "\t", trapezoidal, "\t", error) if abs(error) <= 0.0001: print(n, "nodes") print("By the table above we need n=54 when s=2 to get error less than 1/100.") trapezoidal = CompositeTrapezoidal(f, 2, 2, 54).n() error = (trapezoidal  integral(f(x), x,2,2)).n() print("integral(f, x, 2, 2) ≈", trapezoidal, "by Composite Trapezoid Rule") print() #s=3  print("For Standard Deviation, s, of 3, we have a=3 and b=3") def CompositeTrapezoidal(f, a, b, n): #compute the size of each subinterval h = (ba)/n #here we add all the terms sum = 1/2*f(a) #add the values for all intermediate nodes for i in [1,2,..,n1] : sum += f(a+i*h) #add the last node sum +=1/2*f(b) #multiply by h to complete the computation sum = sum*h #return the sum return sum #define the function we will be using f(x) =(e^((x^2)/2))/(sqrt(2*pi)) #print the top part of the table print ("n", "\t","trapezoidal" , "\t\t", "error") print ("") #call the function several times for power in [21,22,..,55]: n = power trapezoidal = CompositeTrapezoidal(f, 3, 3, n).n() error = (trapezoidal  integral(f(x), x,3,3)).n() if abs(error) >= 0.0001: print (n, "\t", trapezoidal, "\t", error) if abs(error) <= 0.0001: print(n, "nodes") print("By the table above we need n=54 when s=3 to get error less than 1/100.") trapezoidal = CompositeTrapezoidal(f, 3, 3, 54).n() error = (trapezoidal  integral(f(x), x,3,3)).n() print("integral(f, x, 3, 3) ≈", trapezoidal, "by Composite Trapezoid Rule")

Factoring a XOR Expression in an Equation
I have an equation like this:
(i1*x) XOR (i2*x) = r
wherei1
,i2
, andr
are known values.What I need is that same equation but expressed as
i3*x = r
wherei3
should incorporate bothi1
andi2
in some wayBy doing some tests with random values for
x
, I ended up finding that doing(i1*x) XOR (i2*x)
is the same as doingx * (i1 XOR i2) + remainder
.I got this from the distributive property which says that
a*x + b*x = x(a+b)
but this doesn't seem to work with the XOR operator.The workaround with the remainder doesn’t work for me because it adds a second unknown to the equation which I can't even predict.
Is there another way to do this or something I missed?

Syntax error: With multiple conditional function
I am trying to model this equation with multiple conditions but struggling with the correct syntax. I am not sure if I can do it with sympy. My function depends on two variables, one is quiet obvious but the other has three parameters. Here is how i coded it.The Equationequation is also attached
def zp(args): args = T,param p = 1/T yp = log(1p) if length(param) ==2: xi = 0; sig = param(2); mu = param(1); else: xi = param(1); sig = param(2); mu = param(3); break if xi == 0: zp = mu  sig*log(yp); else: zp = mu  (sig/xi)*(1(yp**xi)); break break