I need to find roots for a generalized state space. That is, I have a discrete grid of dimensions grid=AxBx(...)xX
, of which I do not know ex ante how many dimensions it has (the solution should be applicable to any grid.size
) .
I want to find the roots (f(z) = 0
) for every state z
inside grid
using the bisection method. Say remainder
contains f(z)
, and I know f'(z) < 0
. Then I need to
- increase
z
ifremainder
> 0 - decrease
z
ifremainder
< 0
Wlog, say the matrix history
of shape (grid.shape, T)
contains the history of earlier values of z
for every point in the grid and I need to increase z
(since remainder
> 0). I will then need to select zAlternative
inside history[z, :]
that is the "smallest of those, that are larger than z
". In pseudo-code, that is:
zAlternative = hist[z,:][hist[z,:] > z].min()
I had asked this earlier. The solution I was given was
b = sort(history[..., :-1], axis=-1)
mask = b > history[..., -1:]
index = argmax(mask, axis=-1)
indices = tuple([arange(j) for j in b.shape[:-1]])
indices = meshgrid(*indices, indexing='ij', sparse=True)
indices.append(index)
indices = tuple(indices)
lowerZ = history[indices]
b = sort(history[..., :-1], axis=-1)
mask = b <= history[..., -1:]
index = argmax(mask, axis=-1)
indices = tuple([arange(j) for j in b.shape[:-1]])
indices = meshgrid(*indices, indexing='ij', sparse=True)
indices.append(index)
indices = tuple(indices)
higherZ = history[indices]
newZ = history[..., -1]
criterion = 0.05
increase = remainder > 0 + criterion
decrease = remainder < 0 - criterion
newZ[increase] = 0.5*(newZ[increase] + higherZ[increase])
newZ[decrease] = 0.5*(newZ[decrease] + lowerZ[decrease])
However, this code ceases to work for me. I feel extremely bad about admitting it, but I never understood the magic that is happening with the indices, therefore I unfortunately need help.
What the code actually does, it to give me the lowest respectively the highest. That is, if I fix on two specific z
values:
history[z1] = array([0.3, 0.2, 0.1])
history[z2] = array([0.1, 0.2, 0.3])
I will get higherZ[z1]
= 0.3
and lowerZ[z2] = 0.1
, that is, the extrema. The correct value for both cases would have been 0.2
. What's going wrong here?
If needed, in order to generate testing data, you can use something along the lines of
history = tile(array([0.1, 0.3, 0.2, 0.15, 0.13])[newaxis,newaxis,:], (10, 20, 1))
remainder = -1*ones((10, 20))
to test the second case.
Expected outcome
I adjusted the history
variable above, to give test cases for both upwards and downwards. Expected outcome would be
lowerZ = 0.1 * ones((10,20))
higherZ = 0.15 * ones((10,20))
Which is, for every point z
in history[z, :], the next highest previous value (higherZ
) and the next smallest previous value (lowerZ
). Since all points z
have exactly the same history ([0.1, 0.3, 0.2, 0.15, 0.13]
), they will all have the same values for lowerZ
and higherZ
. Of course, in general, the histories for each z
will be different and hence the two matrices will contain potentially different values on every grid point.
z-shaped arrays for the next highest and lowest observation in
history
relative to the current observation, given the current observation ishistory[...,-1:]
This constructs the higher and lower arrays by manipulating the strides of
history
to make it easier to iterate over the observations of each element of z. This is accomplished usingnumpy.lib.stride_tricks.as_strided
and an n-dim generalzed function found at Efficient Overlapping Windows with Numpy - I will include it's source at the endThere is a single python loop that has 200 iterations for
history.shape
of (10,20,x).sliding_window
from Efficient Overlapping Windows with NumpyI compared what you posted here to the solution for your previous post and noticed some differences.
For the smaller z, you said
They said:
For the larger z, you said
They said:
Using the solution for your previous post, I get:
which seems to work