function x=nninit(Phi,Y,Nh)
% Initialize the neural network.

N=size(Y,1);     % Number of data points.
Ni=size(Phi,2);  % Number of input neurons.

% Scale the output data.
ymax=max(Y);
ymin=min(Y);
Ys=2/(ymax-ymin)*(Y-ones(N,1)*ymin)-1;

% Partition the data and estimate weights and biases of the hidden 
% layer.
Nb=floor(N/Nh);   % Number of samples in each partition
Wh=zeros(Ni,Nh);  % Weights
Bh=zeros(1,Nh);   % Biases
for i=1:Nh,
   WB=pinv([Phi(Nb*(i-1)+1:Nb*i,:),ones(Nb,1)])*Ys(Nb*(i-1)+1:Nb*i,:);
   Wh(:,i)=WB(1:Ni,:);
   Bh(:,i)=WB(Ni+1,:);
end;

x(1:Ni*Nh)=Wh(:);                      % Weights of hidden layer
x(Ni*Nh+1:Ni*Nh+Nh)=Bh(:);             % Biases of hidden layer
x(Ni*Nh+Nh+1:Ni*Nh+2*Nh)=randn(Nh,1);  % Weights of output layer
x(Ni*Nh+2*Nh+1)=mean(Y);               % Bias of output layer