[Scilab-users] vector operation replacing for loop

classic Classic list List threaded Threaded
11 messages Options
Heinz Nabielek-3 Heinz Nabielek-3
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] vector operation replacing for loop

Sorry, Scilab friends, I am still not fluid with vector operations.
Can someone rewrite the for loop for me into something much more efficient?

n=1000;
Z=grand(1,n,'nor',0,1);
r=0.9;
V=Z;
for i=2:n;
        V(i)=r*V(i-1)+sqrt(1-r^2)*Z(i);
end;

The transformation generates an autocorrelated (here rho=0.9) normal distribution V from an uncorrelated normal distribution Z and eventually I will need it for very much larger n values....

Heinz
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Guylaine Collewet Guylaine Collewet
Reply | Threaded
Open this post in threaded view
|

Re: vector operation replacing for loop

Hello,

Could you try this?  I think this gives the same result as with the loop
I think that the matrix mat could be probably built in a smarter way

V=zeros(n,1);
V(1)=Z(1);
expo=1:n-1;
vect=[zeros(1,n-1) r.^(expo-1)];
A=makematrix_toeplitz(vect);
mat=A(n:$,1:n-1);
V(2:$)=(Z(1)*r.^expo)'+sqrt(1-r^2)*mat*Z(2:$)';


Guylaine



-----Message d'origine-----
De : users [mailto:[hidden email]] De la part de Heinz
Nabielek
Envoyé : jeudi 1 novembre 2018 09:27
À : Users mailing list for Scilab
Objet : [Scilab-users] vector operation replacing for loop

Sorry, Scilab friends, I am still not fluid with vector operations.
Can someone rewrite the for loop for me into something much more
efficient?

n=1000;
Z=grand(1,n,'nor',0,1);
r=0.9;
V=Z;
for i=2:n;
        V(i)=r*V(i-1)+sqrt(1-r^2)*Z(i);
end;

The transformation generates an autocorrelated (here rho=0.9) normal
distribution V from an uncorrelated normal distribution Z and eventually I
will need it for very much larger n values....

Heinz
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Samuel GOUGEON Samuel GOUGEON
Reply | Threaded
Open this post in threaded view
|

Re: vector operation replacing for loop

In reply to this post by Heinz Nabielek-3
Hello,

Le 01/11/2018 à 09:27, Heinz Nabielek a écrit :
Sorry, Scilab friends, I am still not fluid with vector operations.
Can someone rewrite the for loop for me into something much more efficient?

n=1000;
Z=grand(1,n,'nor',0,1);
r=0.9;
V=Z;
for i=2:n;
	V(i)=r*V(i-1)+sqrt(1-r^2)*Z(i);
end;

The transformation generates an autocorrelated (here rho=0.9) normal distribution V from an uncorrelated normal distribution Z and eventually I will need it for very much larger n values....

You may use filter(), with a feedback component (since V(i) depends on the previous state V(i-1) computed at the previous step).
However, as shown below, a quick trial shows an initial discrepancy between filter() result and yours with the explicit loop.
I don't know why. May be the setting for the initial condition should be carefully considered/tuned...
n=1000;
Z=grand(1,n,'nor',0,1);
r=0.9;
V=Z;
for i=2:n;
    V(i)=r*V(i-1)+sqrt(1-r^2)*Z(i);
end;

y = filter(sqrt(1-r^2),[1 -r], Z);

clf
subplot(3,1,1), plot(Z), ylabel('input Z')
subplot(3,1,2), plot(V), ylabel('V')
d = abs(y-V); d(d==0) = %nan;
subplot(3,1,3), plot2d("nl",d), title('filter() - V')



_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Samuel GOUGEON Samuel GOUGEON
Reply | Threaded
Open this post in threaded view
|

Re: vector operation replacing for loop

Le 01/11/2018 à 13:40, Samuel Gougeon a écrit :
Hello,

Le 01/11/2018 à 09:27, Heinz Nabielek a écrit :
Sorry, Scilab friends, I am still not fluid with vector operations.
Can someone rewrite the for loop for me into something much more efficient?

n=1000;
Z=grand(1,n,'nor',0,1);
r=0.9;
V=Z;
for i=2:n;
	V(i)=r*V(i-1)+sqrt(1-r^2)*Z(i);
end;

The transformation generates an autocorrelated (here rho=0.9) normal distribution V from an uncorrelated normal distribution Z and eventually I will need it for very much larger n values....

You may use filter(), with a feedback component (since V(i) depends on the previous state V(i-1) computed at the previous step).
However, as shown below, a quick trial shows an initial discrepancy between filter() result and yours with the explicit loop.
I don't know why. May be the setting for the initial condition should be carefully considered/tuned...

This is certainly the reason. For i=1, V(i-1) is unknown and must be somewhat provided. The last filter() option zi is not really documented (i have no time to go to the reference to get clear about how "the initial condition relative to a "direct form II transposed" state space representation." works. Trying to provide a value V(0) such that V(1)==Z(1) decreases the initial discrepancy by a factor 10. But it is still non zero.
y = filter(sqrt(1-r^2), [1 -r], Z,  Z(1)*(1-sqrt(1-r^2))/r);


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Samuel GOUGEON Samuel GOUGEON
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] filter(): How to use the last zi argument <= Re: vector operation replacing for loop

Hello,

Is there anyone here that uses to use the last filter() optional argument?

In the example of the thread reminded below:
V = Z;
for i=2:n
   V(i) = a*V(i-1) + b*Z(i);
end

the input is Z, the output is V, a and b are fixed parameters.
The filter() description tells:

Syntax
------
 [y,zf] = filter(B, A, x [,zi])

Parameters
----------
  B  : real vector : the coefficients of the filter numerator in decreasing power order, or
       a polynomial.
  A  : real vector : the coefficients of the filter denominator in decreasing power order,
       or a polynomial.
  x  : real row vector : the input signal
  zi : real row vector of length max(length(a),length(b))-1: the initial condition
       relative to a "direct form II transposed" state space representation. The default
       value is a vector filled with zeros.

So here i assume that zi is V(0) (so actually Z(1)). Then, how to get
[Z(1) V(2:n)]
as from the above loop? I have tried the following, without success:
y = filter(b, [1 -a], Z);
or
y = [Z(1) filter(b, [1 -a], Z(2:$), Z(1))];
or
y
= filter(b, [1 -a], Z, Z(1)*(1-b)/a);

The zi option has no example, and unless filter() is bugged, the way zi works/is taken into account is unclear.

Any hint (before analyzing its hard code..)?

Regards
Samuel


Le 01/11/2018 à 15:16, Samuel Gougeon a écrit :
Le 01/11/2018 à 13:40, Samuel Gougeon a écrit :
Hello,

Le 01/11/2018 à 09:27, Heinz Nabielek a écrit :
Sorry, Scilab friends, I am still not fluid with vector operations.
Can someone rewrite the for loop for me into something much more efficient?

n=1000;
Z=grand(1,n,'nor',0,1);
r=0.9;
V=Z;
for i=2:n;
	V(i)=r*V(i-1)+sqrt(1-r^2)*Z(i);
end;

The transformation generates an autocorrelated (here rho=0.9) normal distribution V from an uncorrelated normal distribution Z and eventually I will need it for very much larger n values....

You may use filter(), with a feedback component (since V(i) depends on the previous state V(i-1) computed at the previous step).
However, as shown below, a quick trial shows an initial discrepancy between filter() result and yours with the explicit loop.
I don't know why. May be the setting for the initial condition should be carefully considered/tuned...

This is certainly the reason. For i=1, V(i-1) is unknown and must be somewhat provided. The last filter() option zi is not really documented (i have no time to go to the reference to get clear about how "the initial condition relative to a "direct form II transposed" state space representation." works. Trying to provide a value V(0) such that V(1)==Z(1) decreases the initial discrepancy by a factor 10. But it is still non zero.
y = filter(sqrt(1-r^2), [1 -r], Z,  Z(1)*(1-sqrt(1-r^2))/r);



_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users



_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
kjubo kjubo
Reply | Threaded
Open this post in threaded view
|

Re: vector operation replacing for loop

In reply to this post by Heinz Nabielek-3


Hello,

precalculate the know values, I gain 40% faster execution time.
Hope helps a bit...

BR

clc, clear, mode(0)

n=1e6
Z=grand(1,n,'nor',0,1);
Z=Z(:);
r=0.9;

tic()
V1=Z;
for ii=2:n;
    V1(ii)=r*V1(ii-1)+sqrt(1-r^2)*Z(ii);
end;
toc()

tic()
V2 = Z;
k1 = sqrt(1-r^2).*Z;
for ii=2:n;
    V2(ii) = r*V2(ii-1) + k1(ii);
end;
toc()

isequal(V1,V2)




--
Sent from: http://mailinglists.scilab.org/Scilab-users-Mailing-Lists-Archives-f2602246.html
_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Heinz Nabielek-3 Heinz Nabielek-3
Reply | Threaded
Open this post in threaded view
|

Re: vector operation replacing for loop

Great. Thanks. Should have thought of it myself...

Still: is there a vector way?
Heinz

> On 07.11.2018, at 21:06, kjubo <[hidden email]> wrote:
>
>
>
> Hello,
>
> precalculate the know values, I gain 40% faster execution time.
> Hope helps a bit...
>
> BR
>
> clc, clear, mode(0)
>
> n=1e6
> Z=grand(1,n,'nor',0,1);
> Z=Z(:);
> r=0.9;
>
> tic()
> V1=Z;
> for ii=2:n;
>    V1(ii)=r*V1(ii-1)+sqrt(1-r^2)*Z(ii);
> end;
> toc()
>
> tic()
> V2 = Z;
> k1 = sqrt(1-r^2).*Z;
> for ii=2:n;
>    V2(ii) = r*V2(ii-1) + k1(ii);
> end;
> toc()
>
> isequal(V1,V2)

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Samuel GOUGEON Samuel GOUGEON
Reply | Threaded
Open this post in threaded view
|

Re: vector operation replacing for loop

Le 08/11/2018 à 20:24, Heinz Nabielek a écrit :
> Great. Thanks. Should have thought of it myself...
>
> Still: is there a vector way?

Didn't you read Guylaine's answer?
She proposes a tricky algorithm. Unfortunately, it cannot be used for
very long Z input, due to memory considerations.
I am afraid that there is no other better "vectorized" solution.

This is why a builtin function like filter() -- that aims to do what you
want -- is required.

BR

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Heinz Nabielek-3 Heinz Nabielek-3
Reply | Threaded
Open this post in threaded view
|

Re: vector operation replacing for loop

Thanks

> On 08.11.2018, at 21:04, Samuel Gougeon <[hidden email]> wrote:
>
> Le 08/11/2018 à 20:24, Heinz Nabielek a écrit :
>> Great. Thanks. Should have thought of it myself...
>>
>> Still: is there a vector way?
>
> Didn't you read Guylaine's answer?
> She proposes a tricky algorithm. Unfortunately, it cannot be used for very long Z input, due to memory considerations.
> I am afraid that there is no other better "vectorized" solution.
>
> This is why a builtin function like filter() -- that aims to do what you want -- is required.
>
> BR
>
> _______________________________________________
> users mailing list
> [hidden email]
> http://lists.scilab.org/mailman/listinfo/users

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Heinz Nabielek-3 Heinz Nabielek-3
Reply | Threaded
Open this post in threaded view
|

Re: filter(): How to use the last zi argument <= Re: vector operation replacing for loop

In reply to this post by Samuel GOUGEON
On 07.11.2018, at 20:57, Samuel Gougeon <[hidden email]> wrote:

>
> Hello,
>
> Is there anyone here that uses to use the last filter() optional argument?
>
> In the example of the thread reminded below:
> V = Z;
> for i=2:n
>    V(i) = a*V(i-1) + b*Z(i);
> end


In my case, a^2 + b^2 =1, if that should make any difference?

Has the problem been sorted out, how to use "filter"?
Heinz





> the input is Z, the output is V, a and b are fixed parameters.
> The filter() description tells:
>
>> Syntax
>> ------
>>  [y,zf] = filter(B, A, x [,zi])
>>
>> Parameters
>> ----------
>>   B  : real vector : the coefficients of the filter numerator in decreasing power order, or
>>        a polynomial.
>>   A  : real vector : the coefficients of the filter denominator in decreasing power order,
>>        or a polynomial.
>>   x  : real row vector : the input signal
>>   zi : real row vector of length max(length(a),length(b))-1: the initial condition
>>        relative to a "direct form II transposed" state space representation. The default
>>        value is a vector filled with zeros.
>
> So here i assume that zi is V(0) (so actually Z(1)). Then, how to get
> [Z(1) V(2:n)]
> as from the above loop? I have tried the following, without success:
> y = filter(b, [1 -a], Z);
> or
> y = [Z(1) filter(b, [1 -a], Z(2:$), Z(1))];
> or
> y = filter(b, [1 -a], Z, Z(1)*(1-b)/a);
>
> The zi option has no example, and unless filter() is bugged, the way zi works/is taken into account is unclear.
>
> Any hint (before analyzing its hard code..)?
>
> Regards
> Samuel
>
>
> Le 01/11/2018 à 15:16, Samuel Gougeon a écrit :
>> Le 01/11/2018 à 13:40, Samuel Gougeon a écrit :
>>> Hello,
>>>
>>> Le 01/11/2018 à 09:27, Heinz Nabielek a écrit :
>>>> Sorry, Scilab friends, I am still not fluid with vector operations.
>>>> Can someone rewrite the for loop for me into something much more efficient?
>>>>
>>>> n=1000;
>>>> Z=grand(1,n,'nor',0,1);
>>>> r=0.9;
>>>> V=Z;
>>>> for i=2:n;
>>>> V(i)=r*V(i-1)+sqrt(1-r^2)*Z(i);
>>>> end;
>>>>
>>>> The transformation generates an autocorrelated (here rho=0.9) normal distribution V from an uncorrelated normal distribution Z and eventually I will need it for very much larger n values....
>>>>
>>>
>>> You may use filter(), with a feedback component (since V(i) depends on the previous state V(i-1) computed at the previous step).
>>> However, as shown below, a quick trial shows an initial discrepancy between filter() result and yours with the explicit loop.
>>> I don't know why. May be the setting for the initial condition should be carefully considered/tuned...
>>
>> This is certainly the reason. For i=1, V(i-1) is unknown and must be somewhat provided. The last filter() option zi is not really documented (i have no time to go to the reference to get clear about how "the initial condition relative to a "direct form II transposed" state space representation." works. Trying to provide a value V(0) such that V(1)==Z(1) decreases the initial discrepancy by a factor 10. But it is still non zero.
>> y = filter(sqrt(1-r^2), [1 -r], Z,  Z(1)*(1-sqrt(1-r^2))/r);
>>
>>
>>
>> _______________________________________________
>> users mailing list
>>
>> [hidden email]
>> http://lists.scilab.org/mailman/listinfo/users
>
>
> _______________________________________________
> users mailing list
> [hidden email]
> http://lists.scilab.org/mailman/listinfo/users

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Samuel GOUGEON Samuel GOUGEON
Reply | Threaded
Open this post in threaded view
|

Re: filter(): How to use the last zi argument <= Re: vector operation replacing for loop

Le 12/11/2018 à 23:57, Heinz Nabielek a écrit :

> On 07.11.2018, at 20:57, Samuel Gougeon <[hidden email]> wrote:
>> Hello,
>>
>> Is there anyone here that uses to use the last filter() optional argument?
>>
>> In the example of the thread reminded below:
>> V = Z;
>> for i=2:n
>>     V(i) = a*V(i-1) + b*Z(i);
>> end
>
> In my case, a^2 + b^2 =1, if that should make any difference?

It does not.

> Has the problem been sorted out, how to use "filter"?

It hasn't been clarified. We likely must go in the code to see what it
does... and document the option.

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users