WO2012008184A1 - 隠れマルコフモデルの推定方法，推定装置および推定プログラム - Google Patents

隠れマルコフモデルの推定方法，推定装置および推定プログラム Download PDF

Info

Publication number
WO2012008184A1
WO2012008184A1 PCT/JP2011/058312 JP2011058312W WO2012008184A1 WO 2012008184 A1 WO2012008184 A1 WO 2012008184A1 JP 2011058312 W JP2011058312 W JP 2011058312W WO 2012008184 A1 WO2012008184 A1 WO 2012008184A1
Authority
WO
WIPO (PCT)
Prior art keywords
probability
subscript
state
value
branch
Prior art date
Application number
PCT/JP2011/058312
Other languages
English (en)
French (fr)
Inventor
泰男 松山
龍之介 林
Original Assignee
学校法人早稲田大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 学校法人早稲田大学 filed Critical 学校法人早稲田大学
Priority to JP2012524465A priority Critical patent/JP5709179B2/ja
Publication of WO2012008184A1 publication Critical patent/WO2012008184A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L15/00—Speech recognition
    • G10L15/08—Speech classification or search
    • G10L15/14—Speech classification or search using statistical models, e.g. Hidden Markov Models [HMMs]
    • G10L15/142—Hidden Markov Models [HMMs]
    • G10L15/144—Training of HMMs
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/29—Graphical models, e.g. Bayesian networks
    • G06F18/295—Markov models or related models, e.g. semi-Markov models; Markov random fields; Networks embedding Markov models

Definitions

  • the present invention relates to a hidden Markov model estimation method, an estimation apparatus, and an estimation program for estimating a hidden Markov model (HMM) that generates a given data series at a high speed.
  • HMM hidden Markov model
  • the algorithm for estimating a hidden Markov model from a given data series was proposed by El. Mr. Baum and L. It is called the Baum-Welch algorithm in the name of L. Welch. Later, it is known that this is the EM (Expectation-Maximization) algorithm specified for the Markovian calculation method.
  • the EM algorithm is a technique for obtaining a probability model that maximizes the likelihood that observed data appears (incomplete data likelihood: incomplete-data likelihood).
  • the log likelihood based on a temporary complete data model ( log-likelihood) from the step of obtaining the conditional expected value based on the observed data (E-step) and the step of updating the model to achieve the maximum of the expected value (M-step) This is accomplished by repeating the two procedures alternately.
  • the EM algorithm is iterative and the use of surrogate functions that are more capable than logarithms is expected to provide a general and fast algorithm structure. Therefore, one of the inventors of the present application has proposed an alpha-EM algorithm including the EM algorithm as a subset in Non-Patent Document 1.
  • This uses alpha-log, which is an extended version of the log (logarithmic) function, and it was expected that there will also be an alpha-EM algorithm version (alpha-HMM re-estimation algorithm) for the HMM re-estimation algorithm. .
  • the present invention provides a hidden Markov model estimation method, estimation apparatus, and estimation program capable of solving the above-described problems and obtaining an HMM estimation algorithm capable of calculating unknown parameters even under conditions other than special conditions.
  • the purpose is to provide.
  • the present invention is a hidden means comprising storage means for storing observation data in time series and estimation means for estimating what kind of probability model the observation data is by calculating unknown parameters of a hidden Markov model.
  • the estimation means includes initial setting means for setting speed-up parameters of the hidden Markov model, and probability of state transition, output, initial state, and likelihood as unknown parameters of the hidden Markov model.
  • Update setting means for updating and setting each quantity and expected value of state transition and output, and each probability quantity and each expected value immediately before update setting by the update setting means, as well as each probability shifted by time before that Using the quantity and each expected value, the observation data read from the storage means, and the acceleration parameter set by the initial setting means And calculating means for calculating each new probability amount and each expected value by applying a micro-approximation, and determining the convergence of the calculation by the calculation means, and if the calculation does not converge, the calculation means calculates A determination means for causing each new probability amount and each expected value to be set and updated by the update setting means and, if the calculation has converged, to output each new probability amount calculated by the calculation means as a final value. It is prepared for.
  • the calculation means calculates the value of the initial state probability ⁇ i
  • ⁇ t + 1 (t + 1 is a subscript of ⁇ ) from the t + 1th state i to the state j is calculated at the tth and t ⁇ 1th times.
  • ⁇ t-1 ij is a subscript of a, t or t-1 is a subscript of ⁇
  • ⁇ t + 1 (t + 1 is a subscript of ⁇ ) at which the state k is output at j is defined as the expected value N bjk
  • the calculation is performed using the formula 65.
  • the update setting means is configured to set the value of the acceleration parameter to 1 ⁇ ⁇ 3.
  • the present invention stores observation data in a time series in a storage means, and what kind of probability model the observation data is, and an estimation means estimates a hidden Markov model by calculating an unknown parameter of the hidden Markov model.
  • the calculation of the unknown parameter includes an initial setting step of setting a speed-up parameter of the hidden Markov model, and probability quantities of state transition, output, initial state, and likelihood as unknown parameters of the hidden Markov model.
  • an update setting step for updating and setting each expected value of state transition and output, and each probability amount and each expected value immediately before the update setting in the update setting step, as well as each probability amount shifted in time before that And using each expected value, setting the observation data read from the storage means and the initial setting step.
  • Each new probability amount and each expected value calculated in the step are set and updated in the update setting step, and if the calculation has converged, each new probability amount calculated in the calculation step is output as a final value.
  • the determination step is performed.
  • the calculation step calculates the value of the initial state probability ⁇ i
  • ⁇ t + 1 (t + 1 is a subscript of ⁇ ) from the + 1st state i to the state j is set as the expected value N aij of the state transition at the tth and t ⁇ 1th times.
  • the calculation step further outputs the expected value of state transition or branching using an arithmetic expression to which l (el) indicating an algorithm repetition index (index) is added. Instead of probability, branch probability, mean value vector, and covariance matrix are calculated.
  • the calculation step adds the index indicating the arrangement number of the observation data and performs the same calculation.
  • the value of the acceleration parameter is set to 1 ⁇ ⁇ 3.
  • the present invention provides a hidden Markov model that causes a computer to function as an estimation means for estimating what kind of probability model the observation data stored in time series in the storage means is by calculating unknown parameters of the hidden Markov model.
  • the estimation means includes initial setting means for setting speed-up parameters of the hidden Markov model, and probability quantities of state transition, output, initial state and likelihood as unknown parameters of the hidden Markov model, and Update setting means for updating and setting each expected value of state transition and output, and not only each probability quantity and each expected value immediately set by the update setting means, but also each probability quantity shifted by time before that While using each expected value, the observation data read from the storage means and the high speed set by the initial setting means And calculating means for calculating each new probability amount and each expected value by applying a small approximation using the parameter, and determining the convergence of the calculation by the calculating means.
  • the calculating means Each new probability amount and each expected value are set and updated by the update setting unit, and if the calculation is converged, the new probability amount by the calculation unit is made to function as a determination unit for determining the final value. Yes.
  • the calculation means calculates the value of the initial state probability ⁇ i
  • ⁇ t + 1 (t + 1 is a subscript of ⁇ ) from the + 1st state i to the state j is set as the expected value N aij of the state transition at the tth and t ⁇ 1th times.
  • the computing means uses the expected value of state transition or branching by an arithmetic expression with l indicating the algorithm iteration index (index), and further substitutes the output probability. Branch probability, mean value vector, and covariance matrix.
  • the calculation means attaches an index that indicates the array of the observation data and performs the same calculation.
  • the update setting means sets the value of the acceleration parameter to 1 ⁇ ⁇ 3.
  • FIG. 6 is a graph showing the correlation between the number of iterations and the likelihood for the convergence speed of the alpha-HMM algorithm.
  • FIG. 6 is a graph showing the correlation between the number of iterations and the likelihood for the convergence speed of the alpha-HMM algorithm.
  • a learning algorithm based on a hidden Markov model is generated starting from optimization of convex divergence.
  • a convex divergence between two probability densities p and q is considered.
  • Y is the K-dimensional Euclidean space.
  • the function f (r) is convex at (0, ⁇ ).
  • the dual function g (r) of the function f (r) needs to satisfy the following equation.
  • Equation 1 the integral in Equation 1 is replaced with the sum. This is the same in all the integrals described below.
  • Equation 2 gives the following equation:
  • o (1) means a high-order infinitesimal term. Note that o (1) is not necessary if a convex function that can be separated by multiplication is used.
  • the following equation is an interesting function, which is a monotone function whose convexity can be adjusted by the parameter c.
  • L (c) (r) is regarded as a regularized logarithm with parameter c. That is, it can be called c-logarithm. However, increasing the parameter c decreases its value. Therefore, here, the ruled logarithm is parameterized by the following relationship.
  • the ⁇ -log is obtained by starting with a convex divergence D f (p ⁇ q) as shown in the following equation.
  • ⁇ divergence is expressed by the following equation.
  • ⁇ ) is the probability density or probability amount of the observation data y in Y parameterized by ⁇ .
  • the set Y is regarded as a set of incomplete data.
  • x ⁇ X (bold) is complete data or augmented data of ideal observation including unknown or lost data.
  • the pdf (probability density function) or pmf (probability mass function) of incomplete data is expressed by the following equation.
  • conditional pdf or pmf is as follows.
  • the alpha-log likelihood ratio of incomplete data in pdf or pmf is as follows:
  • the alpha-EM algorithm is positioned as the starting point of the alpha-HMM re-estimation algorithm.
  • Initial setting Set the parameter ⁇ for the first cycle.
  • E-step The above formula 25 is calculated as a function of ⁇ when ⁇ is given.
  • Update parameters are calculated according to the following formula.
  • Initial setting Set the parameter ⁇ for the first cycle.
  • E-step The above formula 25 is calculated as a function of ⁇ when ⁇ is given.
  • Update parameters are calculated according to the following formula.
  • HMM Hidden Markov Model
  • s i ⁇ is an alphabetic symbol indicating one state in the sum ⁇ within a finite range.
  • the symbol ⁇ used here includes a state transition probability, an output probability, and an initial state probability.
  • Equation 28 the value of the variable string of random observation data shown in Equation 28 is expressed by the following equation.
  • t is an index representing the order of data.
  • Equation 28 The observation data represented by Equation 28 is given as the incomplete data Y described above.
  • Equation B The state transition shown in the following equation is the lost data Z.
  • Equation 38 Let ⁇ and ⁇ be each parameter set of two Markov models. Then, for ⁇ ⁇ 1, the following equation of Equation 38 has the meaning of Equation 39.
  • Equation 39 is established when the inequality of Equation 38 is established.
  • Equation 24 which is a basic equation, the alpha-log likelihood ratio of incomplete data shown in the following equation is set.
  • the characteristic of the equation 41 means that the Markov model ⁇ estimated in the form of the equation 40 generates a sequence y (bold) with a higher probability than the Markov model ⁇ . That is, the iteration of the alpha-EM algorithm in the Markov model class shown in the power equation 33 generates the best hidden Markov model at the cost of inevitable local optimization. Therefore, repeating the equation 38 with the parameters of the Markov model shown in equation 33 is the theoretical version of the HMM algorithm.
  • the main problem with the alpha-EM algorithm is how to provide specific computer algorithms that can be implemented as software.
  • the alpha-HMM of number 38 is a theoretical form, which is a general HMM as well as the log-EM algorithm.
  • Two versions are presented here as alpha-HMMs that can be implemented in software, but the first version is accurate but requires computation of future information. That is, the obtained algorithm is non-causal.
  • the second version particularly proposed in the present invention can be calculated without including future information by using time shift and probability approximation.
  • the integral is the sum.
  • t is a repetition index (index) for alpha-HMM.
  • the parameter set is as shown in the following equation. In the specification, lowercase symbols such as ⁇ , a, and b are shown in capital letters when expressed as a set.
  • N ij (s) (s is bold) is the number of state transitions that transition from i to j in the state s.
  • Formula 50 also includes ⁇ t + 1 on the right side. This will also be resolved in the next version.
  • N bjk (s) (s is bold and jk is a subscript of b) is the number of events with respect to the output probability b jk caused by the sequence s (bold).
  • Formula 53 also includes ⁇ t + 1 on the right side. This will be resolved in the next version.
  • Equation 47 The time shift means that the probability environment for the expected value of the likelihood ratio is reversed. Therefore, the original parameter ⁇ is converted to ⁇ causal .
  • the relationship between ⁇ and ⁇ causal is obtained by the following induction.
  • the following relationship holds in the case of Equation 47.
  • Equation 56 there is a problem in that the right side of the first equation must calculate a function (term) that depends on the same time series ⁇ t + 1 in order to calculate the initial probability ⁇ i of ⁇ t + 1 as the left side. , It included self-adherence that it was impossible to calculate without any change. Therefore, in the second equation, the denominator of the first equation is first divided by the ⁇ power of the term in ⁇ .
  • the concept of time shift is introduced while making a small approximation by Taylor expansion, thereby solving the self-consistency contradiction.
  • the right side does not use a function (term) that depends on ⁇ t + 1 of the simultaneous sequence with the left side. This detailed concept is shown below.
  • Equation 55 and Equation 62 provides a causal and computationally efficient alpha-HMM algorithm that uses past information. It should be noted that this corresponds to the 27 alpha-GEM algorithm.
  • Equation 50 gives the following update equation.
  • ⁇ t (ij is a subscript of a and t is a subscript of ⁇ ) is N ij (s and y are bold) under a conditional probability P (s
  • s) is an expected value of s (s is bold), and is estimated by the count number of state transitions caused by incomplete data y (y is bold).
  • Equation 53 gives the following update equation.
  • ⁇ t (jk is a subscript of b and t is a subscript of ⁇ ) is N bjK (under the conditional probability P (s
  • s) is an expected value of s (s is bold) and is estimated by the count number of state transitions caused by incomplete data y (y is bold).
  • observation data when the observation data is discrete, the case where the observation data has a single array (one observation data) and the case where the observation data has a plurality of arrays (M observation data) will be described. .
  • Examples 3 and 4 when observation data is continuous, a case where observation data has a single array (single observation data) and a case where a plurality of arrays (M observation data) have will be described. .
  • Examples 5 and 6 when the observation data is semi-continuous, the case where the observation data has a single array (one observation data) and the case where the observation data has a plurality of arrays (M observation data) are described. To do.
  • observation data is discrete and continuous, that is, a case where a discrete sequence and a continuous sequence are mixed, and a case where a single sequence is provided and a case where a plurality of sequences are provided will be described together.
  • FIG. 1 shows a processing procedure of a program that enables execution of the above-described alpha-HMM re-estimation algorithm in the first embodiment of the present invention.
  • reference numeral 1 denotes a register as storage means for storing a total of T observation data, and each observation data is stored in time series in order of time ⁇ from 1.
  • y ⁇ represents individual data values
  • y represents a set of T data values
  • observation data with one array is stored in the register 1.
  • the estimation apparatus 10 calculates the unknown parameters of the HMM based on what kind of probability structure (model) each row of data values stored in the register 1 has.
  • Step S1 is a part for setting the value ⁇ of the acceleration parameter described above.
  • Step S2 is a part for determining an initial probability and a convergence determination value. This is performed before repeating steps S3 to S8 described later.
  • the probability structure of the HMM to be finally calculated by the estimation device 10 is expressed by the following equation. This is equal to the aforementioned expression 37.
  • ⁇ S0 (0 is the subscript of s) is the probability in the first state s 0
  • a S ⁇ -1S ⁇ ( ⁇ 1 and ⁇ are subscripts of s) is the state s ⁇ ⁇ 1 at time ⁇ 1.
  • ⁇ 1 is the probability of transition to state s ⁇ at time ⁇
  • b S ⁇ (y ⁇ ) ( ⁇ is a subscript of s) is the probability that y ⁇ is output when transition to state s ⁇ at time ⁇ .
  • the estimation device 10 reads the observation data stored in the register 1 and estimates the initial state probability ⁇ , the state transition probability a, and the output probability b that the observation data is most likely to appear. Then, as shown in the following equation, initial values of the probabilities ⁇ , a, and b are determined.
  • ⁇ 0 represents the combination of initial probabilities ⁇ , a, and b in the 0th iteration as shown in Equation 43, and the estimation device 10 has a state transition probability of transitioning from state i to state j under that condition.
  • ⁇ 0 (0 is a subscript of ⁇ )
  • ⁇ 0 (0 is the subscript of ⁇ ) in which the state k is output in the state j, and the initial state that first becomes the state i
  • ⁇ 0 (0 is a subscript of ⁇ ) is determined.
  • step S2 a convergence determination value for ending the repetition of steps S3 to S9 is determined by determining a convergence range based on the log likelihood. This convergence determination value is used in step S8, and is specifically represented by the following formula 83.
  • Step S3 uses the initial probability first obtained in step S2, and thereafter uses the probability amount and count value updated in step S9, and uses each probability value shown in equation 68 and each equation shown in equation 69. This is the part where the count value is actually set in the memory (not shown) of the estimation apparatus 10.
  • the next step S4 is a part for reading the observation data from the register 1 and calculating the forward probabilities shown in the equations 70 and 71 from the probability amounts of the above equation set in the step S3.
  • ⁇ shown here is a probability value, which is different from the parameter value ⁇ described above.
  • the probability value ⁇ shown in Equation 70 is calculated, and the probability value (likelihood) P (y
  • the next step S5 is a part for reading the observation data from the register 1 and calculating the forward probability shown in the equation 72 from the probability amount of the above equation set in step S3.
  • ⁇ shown here is a probability value and is different from the parameter value ⁇ described above.
  • the probability value ⁇ shown in Equation 72 is calculated.
  • the observation data is read from the register 1 using the probability value ⁇ calculated in step S4 and the probability value ⁇ calculated in step S5, and the count value of the state transition shown in equation 73 and equation 74 This is a part for calculating the output count values shown in FIG.
  • the next step S7 is the condition of the initial state probability ⁇ i
  • step S7 the parameter value ⁇ set in step S1 is used, which is a main feature of the present algorithm.
  • the probability of the subscript t + 1 the self-stickiness could be removed using the time-shifted subscript t and subscript t-1 values.
  • the calculation method is such that the previous past value can be used, the use of the previous past value has led to high speed, and the weight using the past value (speed parameter value) is 1 ⁇ ⁇ ⁇
  • the effect is that it has been expanded to 3.
  • the observation data with one array is stored in the register 1, but the update formula for probability in step S7 in that case is expressed as the following equation.
  • Step S8 is a part for determining convergence of calculation.
  • the likelihood P based on the new probability amount calculated in the immediately preceding step S7 is used, and the calculation in steps S3 to S9 has converged by comparing the likelihood P with the convergence determination value in step S2. It is determined whether or not. Approximately, the number of repetitions may be specified, and when the specified number is reached, it may be determined that the calculation has converged.
  • Step S9 is a part in which the calculated probability amount and event count value are updated for repeated calculation. In this embodiment, two past values are shifted as shown in the following equation.
  • step S3 and step S9 may be combined and the updated probability amount and count value may be set and stored in the memory of the estimation apparatus 10 for the next iteration as they are.
  • step S8 determines whether the calculation has converged. If it is determined in step S8 that the calculation has converged, the process proceeds to step S10, and an HMM using the calculated three sets of probability amounts shown in the following equation is adopted.
  • the estimation device 10 can output the probability structure of the HMM shown in Equation 66 using each value calculated in Equation 79 above.
  • FIG. 2 shows a configuration of the estimation apparatus 10 that realizes the algorithm shown in FIG.
  • the estimation device 10 reads out the observation data stored in the register 1 and the input means 12 incorporating the register 1 as appropriate, and executes the procedure of steps S1 to S10 to obtain the probability structure of the HMM shown in equation 66.
  • the input means 12 of this embodiment has a structure that can store observation data having one array in the register 1, but the estimation apparatus 100 that handles observation data having a plurality of arrays is described in the following embodiment. This will be explained in detail.
  • the hardware configuration of the estimation device 10 can be realized by, for example, a computer including a CPU as an arithmetic processing unit.
  • a program for executing the processing procedure of the estimation apparatus 10 shown in FIG. 1 is stored in a recording medium such as a memory.
  • the program is stored is not limited.
  • the program may be downloaded to a computer via a communication unit.
  • the estimation means 14 includes the speed setting parameter value ⁇ in step S1, the initial setting means 22 for setting the initial probability value and the convergence determination condition in step S2, the setting of the probability amount and the count value in step S3, The update setting means 24 for updating the probability amount and count value in step S9, the forward probability calculation in step S4, the backward probability calculation in step S5, and the state transition count value and output in step S6 from there. And the calculation means 26 for calculating a new probability amount in step S7 and the convergence determination of the calculation in step S8. If the calculation has not converged, the probability set by the update setting means 24 is calculated. While letting you update the quantity and count value If calculations if converged, and includes in step S10, a new probability quantity calculating means 26 is calculated immediately before, and a determination unit 28 for confirming the final parameter values of HMM.
  • the estimation device 10 In addition to the memory 30 for storing the accelerating parameter value, the initial probability value and the convergence determination condition, and the probability amount and count value to be updated and set in a readable / writable manner, the estimation device 10 is necessary. Accordingly, an operation means 32 such as a keyboard and a mouse for operating and inputting a speed-up parameter value and a convergence determination condition is connected to the estimation means 14. The procedure of steps S1 to S10 described above is executed by each of these device configurations.
  • Equation 33 the numerical sequence shown in Equation 33 is used as input data to be stored in the register 1, and its parameters are specified as follows.
  • the initial state was fixed for the following reasons.
  • FIG. 3 shows a convergence tendency with respect to log-likelihood, that is, likelihood P (y
  • the experiment was performed with various values satisfying Equation 60 and Equation 61.
  • step S8 was measured according to the following criteria.
  • P new is the value of the probability P obtained in the current calculation
  • P old is the value of the probability P obtained in the previous calculation.
  • Equation 63, Equation 64, and Equation 65 are exact forms of the log-HMM algorithm. In this sense, the update formulas of Equation 63, Equation 64, and Equation 65 can be regarded as an extended version of the log-HMM algorithm.
  • this embodiment showed a faster HMM algorithm than in the past.
  • This estimation method reflects the underlying alpha-EM algorithm and can be called an improved alpha-HMM algorithm.
  • the improved alpha-HMM algorithm is superior to the conventional alpha-HMM algorithm or log-HMM algorithm. Since the increase in computational complexity is very small, the CPU speed that is directly expected to reduce the number of iterations was realized.
  • ⁇ ⁇ t + 1 and ⁇ ⁇ t ( ⁇ is marked above ⁇ ) indicate that the statistic is estimated from a smaller number of samples (eg, 200 samples are extracted from 1000 samples). . That is, a part of the sample is used to think about the future.
  • This method requires a larger overhead than the update formulas of Equations 63, 64, and 65, but the benefits of speeding up outweigh this extra burden.
  • This method can be used if the user has a priori knowledge that the input sources are well mixed or ergodic. In the experiment, the speedup was not worse than the methods of Equations 63, 64 and 65.
  • the register 1 as a storage means for storing the input observation data in time series and the probability model of the observation data are calculated by calculating the unknown parameters of the HMM.
  • the estimation means 14 includes an initial setting means 22 for setting a value ⁇ of the HMM acceleration parameter, and state transition probabilities a, output probability b, initial state probability ⁇ and likelihood P
  • ⁇ ) probability amounts, state transition expected value N a and output when calculating the expected value N b, by calculating each new probability weights and the expected value by applying a time shift and fine approximation, it is possible to avoid a conflict of conventional self consistent, special conditions Even if it is not below, it is possible to obtain an HMM estimation algorithm capable of calculating unknown parameters. Further, since the format only uses each probability amount and each expected value shifted in time as the accumulated past information, the unknown parameter can be obtained very quickly without taking time for the arithmetic processing.
  • the calculation means 26 here calculates the value of the initial state probability ⁇ i
  • ⁇ t + 1 (t + 1 is a subscript of ⁇ ) from the t + 1th state i to the state j, and the expected value N of the state transition at the tth and t ⁇ 1th times.
  • ⁇ t-1 ij is a subscript of a, t or t-1 is a subscript of ⁇
  • the calculation is performed by the formula 64, and the state k is output in the state j.
  • the update setting means 24 here sets the value of the acceleration parameter to 1 ⁇ ⁇ 3. By doing so, it is possible to reliably increase the speed of calculation of unknown parameters as compared to the conventional method by setting ⁇ to be larger than 1 in the range where ⁇ is less than 3 where repetitive calculations do not diverge.
  • the operational effects of the present embodiment described above are the initial setting step for executing the operation as the initial setting means 22, the update setting step for executing the operation as the update setting means 24, and the operation as the calculating means 26.
  • the HMM estimation method including the calculation step and the determination step for executing the operation as the determination unit 28 is also performed in exactly the same manner, and the HMM that causes the computer to function as the estimation unit 14 Even if it is an estimation program of, it is demonstrated similarly.
  • FIG. 4 shows the processing procedure of the program in the second embodiment of the present invention.
  • reference numeral 101 denotes a register as a storage device for storing a total of T observation data, and each observation data is stored in time series in order of time ⁇ from 1.
  • y ⁇ represents an individual data value
  • y represents a set of T data values
  • n represents an index indicating the number of arrays
  • M ( Integer greater than or equal to 2)
  • This observation data is stored in the register 101.
  • the estimation apparatus 100 estimates and analyzes what kind of probability structure (model) the data values of the M columns stored in the register 101 have according to the following steps S11 to S20. is there.
  • Step S12 is a part for determining an initial probability and a convergence determination value. This corresponds to step S2, and is performed before repeating steps S13 to S18 described later.
  • the estimation apparatus 100 reads the observation data stored in the register 101 and estimates the initial state probability ⁇ , the state transition probability a, and the output probability b at which the observation data is most likely to appear. Then, as shown in the following equation, initial values of the probabilities ⁇ , a, and b are determined.
  • ⁇ 0 represents the combination of initial probabilities ⁇ , a, and b in the 0th iteration, as shown in Equation 43, and the estimating apparatus 100 has a state transition probability of transitioning from state i to state j under that condition.
  • ⁇ 0 (0 is a subscript of ⁇ )
  • ⁇ 0 (0 is the subscript of ⁇ ) in which the state k is output in the state j, and the initial state that first becomes the state i
  • ⁇ 0 (0 is a subscript of ⁇ ) is determined.
  • step S12 a convergence determination value for ending the repetition of steps S13 to S19 is determined by determining a convergence range based on the log likelihood. This convergence determination value is used in step S18, and is specifically expressed by the above equation 83.
  • Step S13 uses the initial probability obtained in Step S12 first, and thereafter uses the probability amount and count value updated in Step S19, and uses each probability value shown in Equation 86 and each equation shown in Equation 87. This is the part that actually sets the count value in the memory (not shown) of the estimation apparatus 100.
  • the next step S14 is a part for reading the observation data from the register 101 and calculating the forward probabilities shown in the equations 88 and 89 from the probability amounts of the above equation set in the step S3.
  • ⁇ shown here is a probability value, which is different from the parameter value ⁇ described above.
  • the probability value ⁇ shown in Equation 88 is calculated, and the probability value (likelihood) P (y (n)
  • the next step S15 is a part for reading the observation data from the register 101 and calculating the forward probability shown in Formula 90 from the probability amount of the above equation set in step S13.
  • ⁇ shown here is a probability value and is different from the parameter value ⁇ described above.
  • the probability value ⁇ shown in Formula 90 is calculated.
  • step S16 using the probability value ⁇ calculated in step S14 and the probability value ⁇ calculated in step S5, the observation data is read out from the register 101, the state transition count value shown in equation 91, and the equation 92 This is a part for calculating the output count values shown in FIG.
  • the next step S17 is a value of the initial state probability ⁇ i
  • step S17 the parameter value ⁇ set in step S11 is used, which is a main characteristic part of the present algorithm.
  • the self-stickiness could be removed using the time-shifted subscript t and subscript t-1 values.
  • the calculation method is such that the previous past value can be used, the use of the previous past value has led to high speed, and the weight using the past value (speed parameter value) is 1 ⁇ ⁇ ⁇
  • the effect is that it has been expanded to 3.
  • M pieces of observation data are stored in the register 101.
  • the probability update formula in step S17 is expressed as the following formula.
  • Step S18 is a part for determining convergence of calculation.
  • the likelihood P based on the new probability amount calculated in the immediately preceding step S17 is used, and the calculation in steps S13 to S19 has converged by comparing the likelihood P with the convergence determination value in step S12. It is determined whether or not. Approximately, the number of repetitions may be specified, and when the specified number is reached, it may be determined that the calculation has converged.
  • Step S19 is a part for updating the calculated probability amount and event count value for repeated calculation. In this embodiment, two past values are shifted as shown in the following equation.
  • the t-th probability amount and count value are updated to the t-1th probability amount and count value, and the t + 1-th probability amount and count value are updated to the t-th probability amount. Updated to probability amount and count value.
  • step S18 determines whether the calculation has converged. If it is determined in step S18 that the calculation has converged, the process proceeds to step S20, and an HMM using the calculated three sets of probability amounts shown in the following equation is adopted.
  • the estimation device 10 can output the probability structure of the HMM shown in Equation 66 using each value calculated in Equation 97 above.
  • FIG. 5 shows the configuration of the estimation apparatus 100 that implements the algorithm shown in FIG.
  • the estimation apparatus 100 reads the observation data stored in the register 101 and the input means 112 incorporating the register 101 as appropriate, and executes the procedure from step S11 to step S20 to obtain the probability structure of the HMM shown in equation 66.
  • the input means 112 of this embodiment has a structure that can store observation data having M arrays in the register 101. Since the hardware configuration of the estimation apparatus 100 is the same as that of the estimation apparatus 10 of the first embodiment, description thereof is omitted here.
  • the estimation unit 114 includes an initial setting unit 122, an update setting unit 124, a calculation unit 126, and a determination unit 128. These correspond to the initial setting means 22, the update setting means 24, the calculation means 26, and the determination means 28 of the first embodiment, respectively, except that the observation data to be handled is arranged in M arrays. Works the same way.
  • the estimation device 110 In addition to the memory 130 for storing the accelerating parameter value, the initial probability value and the convergence determination condition, the updated probability amount and the count value in a readable / writable manner, the estimation device 110 is necessary. Accordingly, an operation means 132 such as a keyboard and a mouse for operating and inputting a speed-up parameter value and a convergence determination condition is connected to the estimation means 114. The procedure of steps S11 to S20 described above is executed by each of these device configurations.
  • the register 101 as a storage means for storing input observation data in time series and the probability model of the observation data are calculated as unknown parameters of the HMM.
  • the estimation means 114 includes an initial setting means 122 for setting the value ⁇ of the HMM acceleration parameter, and a state transition probability a as an unknown parameter of the HMM.
  • Output probability b initial state probability ⁇ and likelihood P (y
  • the observation data read out from the register 101 and the speed-up parameter value ⁇ set by the initial setting means 122 are used to apply a small approximation by Taylor expansion to obtain each new probability amount and each expected value.
  • the calculation means 126 to calculate (refer to Expression 63, Expression 64, Expression 65, Expression 93, Expression 94, Expression 95) and the convergence of the calculation by the operation means 126 are determined.
  • Each of the new probability amounts and each expected value calculated in (5) is set and updated by the update setting means, and if the calculation has converged, the new probability amounts calculated by the calculation means 126 are output as final values. Means 128.
  • ⁇ ) probability amounts, state transition expected value N a and output when calculating the expected value N b, by calculating each new probability weights and the expected value by applying a time shift and fine approximation, it is possible to avoid a conflict of conventional self consistent, special conditions Even if it is not below, it is possible to obtain an HMM estimation algorithm capable of calculating unknown parameters. Further, since the format only uses each probability amount and each expected value shifted in time as the accumulated past information, the unknown parameter can be obtained very quickly without taking time for the arithmetic processing.
  • the computing means 126 here calculates the value of the initial state probability ⁇ i
  • ⁇ t + 1 (t + 1 is a subscript of ⁇ ) from the t + 1th state i to the state j, and the expected value N of the state transition at the tth and t ⁇ 1th times.
  • the update setting means 124 here sets the value of the acceleration parameter to 1 ⁇ ⁇ 3. By doing so, it is possible to reliably increase the speed of calculation of unknown parameters as compared to the conventional method by setting ⁇ to be larger than 1 in the range where ⁇ is less than 3 where repetitive calculations do not diverge.
  • the operational effects of the present embodiment described above execute an initial setting step for executing the operation as the initial setting means 122, an update setting step for executing the operation as the update setting means 124, and an operation as the calculation means 126.
  • the HMM estimation method including the calculation step and the determination step for executing the operation as the determination unit 128 is also performed in exactly the same manner, and the HMM that causes the computer to function as the estimation unit 114 Even if it is an estimation program of, it is demonstrated similarly.
  • observation data is discrete
  • observation data is continuous in a single array
  • the case where the observation data is continuous is a case where each value of the observation data string fluctuates like a wave.
  • the output probabilities b ik and b jk in the model diagram when the observation data in FIG. 6 are discrete are branches that are probability density functions as in the model diagram when the observation data in FIG. 7 is continuous.
  • the probabilities c ik and c jk are replaced.
  • a source data string such as the following expression is given.
  • t is an index representing the order of data.
  • Each y t is a scalar or vector.
  • the above equation is for a single array, but if given by M multiple arrays, the source data sequence is:
  • the symbol ⁇ indicates a set.
  • S S is bold
  • S means a state transition sequence as shown below.
  • Each probability is as shown in the following formulas 102 to 105.
  • Equation 102 to 105 The set of probabilities shown in Equations 102 to 105 above is collectively expressed as ⁇ in the following equation.
  • the probabilistic data structure is interpreted as follows.
  • the incomplete data is y (y is bold)
  • the estimated disappearance data is s (s is bold)
  • the complete data has several hundred probability quantities and several 106 parameters, It is defined as x (x is bold) shown in the following formula.
  • Alpha-EM algorithm HMM interpretation with incomplete data, missing data, and complete data is comparable to the EM algorithm.
  • the purpose of this paper is to find a new HMM estimation algorithm, so it shows the path starting from the alpha-EM algorithm.
  • observation data y (y is bold) is parameterized by ⁇ representing the entire parameter to be estimated
  • ⁇ ) (y is bold) is a probability density or a probability amount.
  • x ⁇ X (x and X are bold) be complete data that is an ideal observation result including lost data, or expanded data.
  • An incomplete data probability density function (pdf) or probability quantity function (pmf) is expressed by the following equation.
  • the alpha-EM algorithm uses the alpha logarithm as shown below.
  • ⁇ and ⁇ indicate the parameters of the old model and the new model for the number 106 in the iteration maximization step. Then, the basic equation of alpha-EM algorithm is obtained as
  • D ( ⁇ ) is two conditional probabilities P x
  • ⁇ ) alpha and divergence between x and y are bold.
  • Q function is important.
  • E indicates the objective function for maximization represented by the quantity.
  • the alpha log likelihood ratio of the incomplete data which is the left side of the expression 114 is also positive in the range of a ⁇ 1. Therefore, the alpha-EM algorithm and its variant alpha-GEM algorithm are shown below.
  • Equation 106 The initial value of Equation 106 is selected and used as ⁇ .
  • E-Step The calculation of Formula 115 is executed.
  • M-step The update parameter is calculated by the following formula.
  • Arg max is a symbol that means a variable that gives the maximum value.
  • the following equation expresses Q ( ⁇ ) (x
  • the approximate version of the alpha-HMM algorithm is the alpha-GEM algorithm.
  • the basic characteristics of the alpha-HMM estimation algorithm due to the advantages of the alpha-EM algorithm having Formula 114 will be described.
  • the update equation of the state transition probability a ij will be described. Since the state transition probability a ij needs to be a probability amount after the update, it is necessary to use Lagrange's undetermined coefficient method (Lagrange multiplier). Therefore, differentiation for calculating the maximum value is performed as in the following equation.
  • Lagrange multiplier Lagrange multiplier
  • N ij (s) (s is bold) is the number and position of state transitions from i to j.
  • N ij (s) (s is bold) is the number and position of state transitions from i to j.
  • N T the order of O
  • ⁇ causal is hereinafter referred to as ⁇ .
  • Equation 125 which is a causal approximation
  • sequence extension of Equation 127 give the following update equations for transition probabilities. Note that the above formulas 122 to 124, which are update formulas of these probabilities, correspond to the formulas 75 to 77 of the first embodiment, respectively.
  • ⁇ l + 1 (l + 1 is a subscript of ⁇ ) of a single-array discrete type alpha-HMM is as follows.
  • ⁇ l + 1 (l + 1 is a subscript of ⁇ ) of a single-array discrete type alpha-HMM is as follows.
  • ⁇ l + 1 (l + 1 is a subscript of ⁇ ) of a single array discrete type alpha-HMM is as follows.
  • each probability amount is calculated based on the following calculation using a finite number of past information.
  • the speed-up parameter ⁇ is a design parameter selected by the user, and has an effect of weighting the previous repeated value to the repeated value based on the current value.
  • the weighting and addition to the past values can be performed retroactively to a finite past. That is, ⁇ ⁇ ⁇ 0 may be set as a design parameter and replaced as in the following equation.
  • Expression 131 and Expression 134 can be expanded as the following expression (for example, Taylor expansion).
  • o (1) is a symbol meaning high-order infinitesimal
  • the first term on the right side of the above equation corresponds to a conventional hidden Markov model estimation algorithm.
  • the second term on the right side corresponding to the momentum term has the meaning of a speed-up term approximated by expansion.
  • the performance curve as shown in FIG. 3 at the rising portion where the calculation converges, wavy up and down fluctuations are likely to occur as the speed increases, but the fluctuations are reduced.
  • S (S is bold) is a set of M state transition sequences as shown in the following equation.
  • the Q function of a plurality of arrays is as follows.
  • P is the probability of the Markov process shown by the following equation.
  • initial state probability ⁇ s0 (n) (0 is a subscript of s, (n) is a subscript), state transition probability a st-1 (n) st ( n) (t-1 and t are subscripts of s, (n) is a subscript of subscripts), and output probability b st (n) (y t (n) ) (t Is a subscript of s, and (n) is a superscript) is independent of the series index (index) n.
  • MLE maximum likelihood estimation
  • the output probability b stkt (y t) (stkt of t is, s respectively, subscript subscripts k) in the branch probability C that identifies the probability of transition to the k t -th branch stkt (t of stkt is, s respectively, subscript subscript k) is the probability density function of y t.
  • the output probability b jk (y t ) at which the state k is output in the state j is expressed by the following equation.
  • ⁇ jk is an average value vector
  • ⁇ jk is a covariance matrix
  • the average value vector indicates the average value of the probability density function of the branching probability (the same applies hereinafter).
  • the output probability density function (pdf) in the state j is expressed by the following equation so that b j (y t ) becomes pdf.
  • N (y t ; ⁇ jk , ⁇ jk ) is a multidimensional normal probability density function (multidimensional Gaussian probability density function) having a mean value vector ⁇ jk and ⁇ jk as a covariance matrix, and y t Is a distribution variable.
  • FIG. 6 illustrates the case of a discrete alphabet.
  • FIG. 7 shows the case of a Gaussian mixture model (bottom), that is, the case of a continuous alphabet. Referring to these drawings, it can be seen that the branch probability c jk arrow in FIG. 7 corresponds to the output probability b jk arrow in FIG.
  • the update equations for the initial probability and the state transition probability are the same as those in the equations 128 and 130, respectively.
  • the update equation of the branch probability c jk is obtained by the Lagrange multiplier coefficient method (Lagrange multiplier) as in the case of the state transition probability a ij .
  • the following equation is an update equation for the mean value vector ⁇ jk .
  • the following non-causal equation is obtained as a direct derivative of Formula 152 for ⁇ jk .
  • Equation 156 shows that the past information of ⁇ l-1 can be used sufficiently.
  • Updating the covariance matrix requires matrix differentiation.
  • the following non-causal equation is obtained by differentiating the number 152 of the Q function with respect to the inverse matrix ⁇ ⁇ 1 jk of the covariance matrix (jk is a subscript of ⁇ ).
  • the update formula of the covariance matrix has a form of past information that can be used effectively and effectively.
  • ⁇ l; 1 represents the covariance matrix of the output when shifting to the branch k in the (l + 1) th state j (hereinafter the same).
  • FIG. 8 shows the processing procedure of the program in the third embodiment of the present invention.
  • the overall flow of this flowchart is substantially the same as in the first embodiment.
  • the present embodiment is different in that the observation data y is single-sequence continuous series data appearing as continuous multivariable observation results.
  • unknown parameters that are set as initial values and repeatedly calculated include initial state probabilities, state transition probabilities, branch probabilities, average value vectors, and covariance matrices, and therefore the calculation formulas in each step are also different.
  • Equation 149 The probability structure of the HMM to be finally calculated by the estimation apparatus 10 is expressed by equation 149 instead of equation 66.
  • the initial value of the unknown parameter determined in step S32 is expressed by the following equation instead of Equation 67.
  • step S32 The method for determining the convergence determination value in step S32 is the same as in equation 83.
  • step S33 the initial value of the unknown parameter first determined in step S32 is used, and thereafter, using the unknown parameter and the count value updated in step S39, each probability value shown in the following equation 163, Each count value shown in Formula 164 is actually set in a memory (not shown) of the estimation apparatus 10.
  • the index (index) is set by adding l of the repetition of the algorithm instead of t in the data order.
  • steps S34 to S36 are performed by replacing a part of the symbols of the formulas 70 to 74 in the first embodiment. Specifically, the output probability b j. y ⁇ + 1 ( ⁇ + 1 is a subscript of y) and branch probability c j. y ⁇ + 1 ( ⁇ + 1 is a subscript of y) and t (index) is replaced with l.
  • step S37 the initial state probability initial state probability ⁇ i
  • the speed-up parameter value ⁇ is used, and in calculating the probability of the portion of the subscript l + 1, self-adherence using the time-shifted subscript l and subscript l-1 values is used.
  • the result is a calculation method that can use the previous past value, the use of the previous past value has led to high speed, and the past value is used.
  • the fact that the weight (value of the speed-up parameter) is expanded to 1 ⁇ ⁇ ⁇ 3 can be cited as an effect as in the first embodiment.
  • the observation data with one array is stored in the register 1, and the update formula for probability in step S37 in that case is as described above. That is, the initial state probability update formula is Formula 130, the state transition probability update formula is Formula 128, the branch probability update formula is Formula 154, the mean value vector update formula is Formula 156, and the covariance matrix update formula is Formula 158. ⁇ 160.
  • step S38 as in the first embodiment, the likelihood P based on the new probability amount calculated in step 37 is used, and the likelihood P is compared with the convergence determination value in step S32. It is determined whether or not the calculation of S39 has converged. If not converged, the process proceeds to step S39, the unknown parameter including the probability amount and the event count value are updated, and the calculation is repeated. At this time, two past values are shifted as shown in the following equation.
  • step S33 and step S39 may be combined and the updated probability amount and count value may be set and stored in the memory of the estimation apparatus 10 for the next iteration as they are.
  • step S38 determines whether the calculation has converged. If it is determined in step S38 that the calculation has converged, the process proceeds to step S40, and the HMM using the calculated five sets of probability amounts shown in the following equation is adopted.
  • the estimation device 10 can output the probability structure of the HMM shown in Equation 149 using each value calculated in Equation 166.
  • the configuration of the estimation device that realizes the algorithm shown in FIG. 8 is the same as that of the estimation device 10 of the first embodiment shown in FIG.
  • the register 1 as a storage means for storing the input observation data in time series and the probability model of the observation data are calculated by calculating the unknown parameters of the HMM.
  • the estimation means 14 includes an initial setting means 22 for setting a value ⁇ of the HMM acceleration parameter, and state transition probabilities a, Each probability amount of branch probability c, mean value vector ⁇ , covariance matrix ⁇ , initial state probability ⁇ , and likelihood P (y
  • Update setting means 24 to be used, and not only the respective probability amounts and expected values immediately before being updated by the update setting means 24 but also the respective probability amounts and expected values shifted in time before that (Equation 16) 164, and the observation data read from the register 1 and the speed-up parameter value ⁇ set by the initial setting means 22 and applying a small approximation by Taylor expansion to each new probability amount and Each expectation value is calculated (see Expression 128, Expression 130, Expression 154, Expression 156, Expression 158 to Expression 160), and the convergence of the calculation by the operation means 26 is determined. Each new probability amount and each expected value calculated by the calculation means 26 are set and updated by the update setting means, and if the calculation has converged, each new probability amount calculated by the calculation means 26 is converted into a final value. And determining means 28 for outputting as follows.
  • the calculation means 26 here calculates the value of the initial state probability ⁇ i
  • ⁇ l + 1 (l + 1 is a subscript of ⁇ ) from the state i + 1 to the state j, and the expected value N of the state transition at the lth and l ⁇ 1th times
  • ⁇ l-1 ij is a subscript of a, l or l-1 is a subscript of ⁇
  • ⁇ l + 1 specifying the probability of transition to the branch transitioning to branch k when the output probability b jk
  • ⁇ l + 1 (l + 1 is a subscript of ⁇
  • the update setting means 24 here sets the value of the acceleration parameter to 1 ⁇ ⁇ 3. By doing so, it is possible to reliably increase the speed of calculation of unknown parameters as compared to the conventional method by setting ⁇ to be larger than 1 in the range where ⁇ is less than 3 where repetitive calculations do not diverge.
  • the operational effects of the present embodiment described above are the initial setting step for executing the operation as the initial setting means 22, the update setting step for executing the operation as the update setting means 24, and the operation as the calculating means 26.
  • the HMM estimation method including the calculation step and the determination step for executing the operation as the determination unit 28 is also performed in exactly the same manner, and the HMM that causes the computer to function as the estimation unit 14 Even if it is an estimation program of, it is demonstrated similarly.
  • this embodiment can be applied to continuous observation data, it can be applied not only to speech recognition by a computer and robot motion recognition, but also to speech synthesis and robot motion generation.
  • the combination of sets of unknown parameters that are repeatedly updated is expressed by Formula 161, and an update equation for alpha-HMM of continuous symbols is obtained.
  • the update equations for the initial state probability and the state transition probability are the same as in equations 142 to 145, respectively.
  • the set of output update equations is different from the discrete symbol case. We need branching probabilities, mean vectors, and covariance matrix update equations.
  • the branch probability can be calculated using Lagrange's undetermined coefficient method as in the case of the initial state probability and the state transition probability. Then, the following update equation is obtained.
  • the covariance matrix update equation is a matrix derivative with respect to ⁇ ⁇ 1 jk
  • G (n) ⁇ l ( ⁇ l is a subscript of G) is the following equation.
  • FIG. 9 shows the processing procedure of the program in the fourth embodiment of the present invention.
  • the overall flow of this flowchart is substantially the same as in the case of the third embodiment.
  • the present embodiment is different in that the observation data y is continuous sequence data of a plurality of arrays that appear as continuous multivariable observation results.
  • step S42 The method for determining the convergence determination value in step S42 is the same as in equation 83.
  • step S43 the initial value of the unknown parameter first determined in step S42 is used, and thereafter, the unknown parameter and the count value updated in step S49 are used.
  • the index (index) is set by adding l of the repetition of the algorithm instead of t in the data order.
  • steps S44 to S46 are performed by replacing a part of the symbols in the expressions 88 to 92 in the second embodiment. Specifically, the output probability b j. y ⁇ + 1 ( ⁇ + 1 is a subscript of y) and branch probability c j. y ⁇ + 1 ( ⁇ + 1 is a subscript of y) and t (index) is replaced with l.
  • step S47 the initial state probability ⁇ i
  • +1 is the subscript of ⁇
  • ⁇ l + 1 (l + 1 is the subscript of ⁇ ) from state i to state j
  • state k is output in state j.
  • the speed-up parameter value ⁇ is used, and in calculating the probability of the portion of the subscript l + 1, self-adherence using the time-shifted subscript l and subscript l-1 values is used.
  • the result is a calculation method that can use the previous past value, the use of the previous past value has led to high speed, and the past value is used.
  • the fact that the weight (value of the speed-up parameter) is expanded to 1 ⁇ ⁇ ⁇ 3 can be cited as an effect as in the second embodiment.
  • M pieces of observation data are stored in the register 101.
  • the probability update formula in step S47 is as described above. That is, the initial state probability update formula is Formula 142, the state transition probability update formula is Formula 143 to Formula 145, the branch probability update formula is Formula 167 to Formula 169, the average value vector update formula is Formula 170, Formula 171; The covariance matrix update formulas are 172 and 173.
  • step S48 in the same manner as in the second embodiment, the likelihood P based on the new probability amount calculated in step S47 is used, and the likelihood P is compared with the convergence determination value in step S42. It is determined whether or not the calculation of S49 has converged. If not converged, the process proceeds to step S49 where the unknown parameter including the probability amount and the event count value are updated and the calculation is repeated. At this time, two past values are shifted as shown in the following equation.
  • the 1st probability amount and count value are updated to the (l-1) th probability amount and count value, and the l + 1st probability amount and count value are updated to the lth time. Updated to probability amount and count value.
  • step S48 determines whether the calculation has converged. If it is determined in step S48 that the calculation has converged, the process proceeds to step S50, and the HMM using the calculated five sets of probability amounts shown in the following equation is adopted.
  • the estimation apparatus 101 can output the probability structure of the HMM shown in Equation 149 using each value calculated in Equation 177.
  • the configuration of the estimation device that realizes the algorithm shown in FIG. 9 is the same as that of the estimation device 101 of the second embodiment shown in FIG.
  • the register 101 as a storage means for storing input observation data in time series and the probability model of the observation data are calculated as unknown parameters of the HMM.
  • the estimation means 114 includes an initial setting means 122 for setting the value ⁇ of the HMM acceleration parameter, and a state transition probability a as an unknown parameter of the HMM.
  • Branch probability c mean value vector ⁇ , covariance matrix ⁇ , initial state probability ⁇ , likelihood P (y
  • Update setting means 124 to be set, not only the respective probability amounts and expected values immediately before being updated and set by the update setting means 124, but also the respective probability amounts and respective time shifts before that A waiting value is used (see Equations 174 and 175), and observation data read from the register 101 and a speed-up parameter value ⁇ set by the initial setting means 122 are used to apply a small approximation by Taylor expansion. Calculation of each new probability amount and each expected value (see Equation 142 to Equation 145, Equation 167 to Equation 173), and the convergence of the calculation by the operator 126 is determined, and the calculation must be converged.
  • each new probability amount and each expected value calculated by the calculation unit 126 are set and updated by the update setting unit, and if the calculation has converged, each new probability amount calculated by the calculation unit 126 is finally updated.
  • Determination means 128 for outputting the value as a value.
  • the computing means 126 here calculates the value of the initial state probability ⁇ i
  • ⁇ l-1 ij is a subscript of a, l or l-1 is a subscript of ⁇
  • calculation is performed using the equations 143 to 145, and the state is j
  • the update setting means 124 here sets the value of the acceleration parameter to 1 ⁇ ⁇ 3. By doing so, it is possible to reliably increase the speed of calculation of unknown parameters as compared to the conventional method by setting ⁇ to be larger than 1 in the range where ⁇ is less than 3 where repetitive calculations do not diverge.
  • the operational effects of the present embodiment described above are the initial setting step for executing the operation as the initial setting means 122, the update setting step for executing the operation as the update setting means 124, and the operation as the calculating means 126.
  • An HMM estimation method including a calculation step and a determination step for executing an operation as the determination unit 128 is also performed in exactly the same manner, and an HMM that causes a computer to function as such an estimation unit 114. Even if it is an estimation program of, it is demonstrated similarly.
  • the present embodiment can be applied to continuous series of observation data, it can be applied not only to speech recognition by a computer and robot motion recognition, but also to speech synthesis and robot motion generation.
  • FIG. 7 the following is realized by reviewing the graph structure of a Gaussian mixture alpha-HMM (also log-HMM).
  • A In a Gaussian mixed HMM, each Gaussian pdf depends on the arriving state j. All N ⁇ K Gaussian density learning requires a variety of long training sequences.
  • the semi-continuous alpha-HMM model is transformed from the equation 149 into the following equation.
  • the update equation of semi-continuous alpha-HMM is as follows. [Initial state probability of single-sequence semi-continuous alpha-HMM] The update formula is Equation 130. [State transition probability of single-sequence semi-continuous alpha-HMM] The update formula is Expression 128. [Branch probability of single-sequence semi-continuous alpha-HMM] The update formula is 154.
  • the processing procedure of the program in the fifth embodiment of the present invention is almost the same as that in the third embodiment, it will be described with reference to FIG.
  • the present embodiment is different in that the observation data y is single-sequence semi-continuous series data appearing as continuous multivariable observation results.
  • step S32 The method for determining the convergence determination value in step S32 is the same as in equation 83.
  • step S33 the initial value of the unknown parameter first determined in step S32 is used, and thereafter, using the unknown parameter and the count value updated in step S39, each probability value shown in the following equation 184, Each count value shown in Expression 185 is actually set in a memory (not shown) of the estimation apparatus 10.
  • steps S34 to S36 are performed by replacing a part of the symbols of the formulas 70 to 74 of the first embodiment and applying the same as in the third embodiment.
  • step S37 the initial state probability ⁇ i
  • the speed-up parameter value ⁇ is used, and in calculating the probability of the portion of the subscript l + 1, self-adherence using the time-shifted subscript l and subscript l-1 values is used.
  • the result is a calculation method that can use the previous past value, the use of the previous past value has led to high speed, and the past value is used.
  • the fact that the weight (value of the speed-up parameter) is expanded to 1 ⁇ ⁇ ⁇ 3 can be cited as an effect as in the third embodiment.
  • the observation data with one array is stored in the register 1, and the update formula for probability in step S37 in that case is as described above. That is, the initial state probability update formula is Formula 130, the state transition probability update formula is Formula 128, the branch probability update formula is Formula 154, the mean value vector update formula is Formula 179, and the covariance matrix update formula is Formula 180. Is a number 182.
  • step S38 as in the third embodiment, the likelihood P based on the new probability amount calculated in step 37 is used, and the likelihood P is compared with the convergence determination value in step S32. It is determined whether or not the calculation of S39 has converged. If not converged, the process proceeds to step S39, the unknown parameter including the probability amount and the event count value are updated, and the calculation is repeated. At this time, two past values are shifted as shown in the following equation.
  • step S33 and step S39 may be combined and the updated probability amount and count value may be set and stored in the memory of the estimation apparatus 10 for the next iteration as they are.
  • step S38 determines whether the calculation has converged. If it is determined in step S38 that the calculation has converged, the process proceeds to step S40, and the HMM using the calculated five sets of probability amounts shown in the following equation is adopted.
  • the estimation device 10 can output the probability structure of the HMM shown in Equation 178 using each value calculated in Equation 187.
  • the configuration of the estimation device that realizes the algorithm shown in FIG. 8 is the same as that of the estimation device 10 of the third embodiment shown in FIG.
  • the unknown parameter of the HMM is calculated based on the register 1 as a storage means for storing the input observation data in time series and what kind of probability model the observation data is.
  • the estimation means 14 includes an initial setting means 22 for setting the value ⁇ of the HMM acceleration parameter, and a state transition probability a , Branch probability c, mean value vector ⁇ , covariance matrix ⁇ , initial state probability ⁇ , likelihood P (y
  • the update setting means 24 to be set and not only the respective probability amounts and expected values immediately before being updated by the update setting means 24 but also the respective probability amounts and expected values shifted in time before that are used ( 184 and Equation 185), and using the observation data read out from the register 1 and the speed-up parameter value ⁇ set by the initial setting means 22 and applying a small approximation by Taylor expansion, each new probability amount And calculating each expected value (see Equation 130, Equation 128, Equation 154, Equation 179 to Equation 182), and the convergence of the calculation by the operator 26 is determined. Each new probability amount and each expected value calculated by the means 26 are set and updated by the update setting means, and if the calculation has converged, each new probability amount calculated by the calculation means 26 is output as a final value. And determination means 28 for making it.
  • the calculation means 26 here calculates the value of the initial state probability ⁇ i
  • the update setting means 24 here sets the value of the acceleration parameter to 1 ⁇ ⁇ 3. By doing so, it is possible to reliably increase the speed of calculation of unknown parameters as compared to the conventional method by setting ⁇ to be larger than 1 in the range where ⁇ is less than 3 where repetitive calculations do not diverge.
  • the operational effects of the present embodiment described above are the initial setting step for executing the operation as the initial setting means 22, the update setting step for executing the operation as the update setting means 24, and the operation as the calculating means 26.
  • the HMM estimation method including the calculation step and the determination step for executing the operation as the determination means 28 is also performed in exactly the same manner, and the HMM that causes the computer to function as such an estimation means 14. Even if it is an estimation program of, it is demonstrated similarly.
  • the present embodiment can be applied to continuous series of observation data, it can be applied not only to speech recognition by a computer and robot motion recognition, but also to speech synthesis and robot motion generation.
  • the update equation for multiple sequences in this case is obtained by limiting the state dependence of the mean vector and covariance.
  • the update equations for the initial state probability, the state transition probability, and the branch probability are the same as those in Equations 142 to 145 and Equations 167 to 169. Therefore, the following equation is obtained.
  • a bar (-) is added above each symbol to distinguish the symbol from the case without a bar (-). This is because, for example, the mathematical expression is different between Expression 171 and the following Expression 189, and the same symbol cannot be used. That is, the bar (-) itself has no special meaning.
  • the covariance matrix update equation is obtained as follows by removing the state dependency.
  • the processing procedure of the program in the sixth embodiment of the present invention is almost the same as that in the fourth embodiment, and will be described with reference to FIG.
  • the present embodiment is different in that the observation data y is a plurality of sequences of semi-continuous series data appearing as continuous multi-variable observation results.
  • step S42 The method for determining the convergence determination value in step S42 is the same as in equation 83.
  • step S43 the initial value of the unknown parameter first determined in step S42 is used, and thereafter, using the unknown parameter and the count value updated in step S49, each probability value shown in the following equation 192, Each count value shown in Equation 193 is actually set in a memory (not shown) of the estimation apparatus 100.
  • steps S44 to S46 are performed by replacing some of the symbols of the formulas 88 to 92 in the second embodiment in the same manner as in the fourth embodiment.
  • step S47 an initial state probability ⁇ i
  • the speed-up parameter value ⁇ is used, and in calculating the probability of the portion of the subscript l + 1, self-adherence using the time-shifted subscript l and subscript l-1 values is used.
  • the result is a calculation method that can use the previous past value, the use of the previous past value has led to high speed, and the past value is used.
  • the fact that the weight (value of the speed-up parameter) is expanded to 1 ⁇ ⁇ ⁇ 3 can be cited as the effect as in the fourth embodiment.
  • M pieces of observation data are stored in the register 101.
  • the probability update formula in step S47 is as described above. That is, the initial state probability update formula is Formula 142, the state transition probability update formula is Formula 143 to Formula 145, the branch probability update formula is Formula 167 to Formula 169, the average value vector update formula is Formula 188, Formula 189, Covariance matrix update formulas are 190 and 191, respectively.
  • step S48 as in the fourth embodiment, the likelihood P based on the new probability amount calculated in step S47 is used, and the likelihood P is compared with the convergence determination value in step S42. It is determined whether or not the calculation of S49 has converged. If not converged, the process proceeds to step S49 where the unknown parameter including the probability amount and the event count value are updated and the calculation is repeated. At this time, two past values are shifted as shown in the following equation.
  • the 1st probability amount and count value are updated to the (l-1) th probability amount and count value, and the l + 1st probability amount and count value are updated to the lth time. Updated to probability amount and count value.
  • step S48 determines whether the calculation has converged. If it is determined in step S48 that the calculation has converged, the process proceeds to step S50, and the HMM using the calculated five sets of probability amounts shown in the following equation is adopted.
  • the estimation apparatus 101 can output the probability structure of the HMM shown in Equation 178 using each value calculated in Equation 195.
  • the configuration of the estimation device that realizes the algorithm shown in FIG. 9 is the same as that of the estimation device 101 of the second embodiment shown in FIG.
  • the register 101 as a storage means for storing input observation data in time series and the probability model of the observation data are calculated as unknown parameters of the HMM.
  • the estimation means 114 includes an initial setting means 122 for setting the value ⁇ of the HMM acceleration parameter, and a state transition probability a , Branch probability c, mean value vector ⁇ ( ⁇ is attached to ⁇ ), covariance matrix ⁇ ( ⁇ is attached to ⁇ ), initial state probability ⁇ , and likelihood P (y
  • Update amount setting means 124 for updating and setting the probability amount, the expected value N a for state transition and the expected value N c for branch, and not only the probability amount and each expected value immediately before the update setting means 124 for update setting, Than Each probability amount and each expected value shifted in the previous time are used (see Equations 192 and 193), and the observation data read from the register 101 and the speed-up parameter value ⁇ set by the initial setting means 122 are used (see Equations 192 and 193), and the observation data read from the register 101 and the speed-up
  • each new probability amount and each expected value by applying a micro-approximation by Taylor expansion (see Equations 142 to 145, Equations 167 to 169, Equations 188 to 191), The convergence of the calculation by the means 126 is determined. If the calculation has not converged, each new probability amount and each expected value calculated by the computing means 126 are set and updated by the update setting means, and if the calculation has converged. And a determination unit 128 that outputs each new probability amount calculated by the calculation unit 126 as a final value.
  • the state transition probability a the branching probability c, the average value vector ⁇ ( ⁇ is attached to ⁇ ), and the covariance matrix ⁇ ( ⁇ is attached to ⁇ )
  • ⁇ ) the state transition expectation value N a
  • the computing means 126 here calculates the value of the initial state probability ⁇ i
  • ⁇ l-1 ij is a subscript of a, l or l-1 is a subscript of ⁇
  • calculation is performed using the equations 143 to 145, and the state is j the output probability b jk
  • the update setting means 124 here sets the value of the acceleration parameter to 1 ⁇ ⁇ 3. By doing so, it is possible to reliably increase the speed of calculation of unknown parameters as compared to the conventional method by setting ⁇ to be larger than 1 in the range where ⁇ is less than 3 where repetitive calculations do not diverge.
  • the operational effects of the present embodiment described above are the initial setting step for executing the operation as the initial setting means 122, the update setting step for executing the operation as the update setting means 124, and the operation as the calculating means 126.
  • An HMM estimation method including a calculation step and a determination step for executing an operation as the determination unit 128 is also performed in exactly the same manner, and an HMM that causes a computer to function as such an estimation unit 114. Even if it is an estimation program of, it is demonstrated similarly.
  • the present embodiment can be applied to continuous series of observation data, it can be applied not only to speech recognition by a computer and robot motion recognition, but also to speech synthesis and robot motion generation.
  • D g is a set whose parts do not overlap.
  • G is the density of the part. If D g corresponds to a particular subclass of a continuous alphabet, this can be considered that the alphabet is accompanied by discrete symbols. In other words, in each section where data is continuous, the sum is taken, and the sum is taken for each sum. Also in the case of this discrete and continuous alphabet, there are cases of a single array and a plurality of arrays as in the above embodiment. Therefore, this application has a total of eight types of alpha-HMM.
  • the discrete value may be considered as a label when a group that outputs continuous values is grouped. That is, the data series is shown as follows:
  • c t g ⁇ .
  • Equation 142 In the case of multiple arrays of discrete continuous mixed series [Initial state probability of alpha-HMM for discrete continuous mixed series of multiple arrays] The initial state probability is the same as Equation 142. [Probability of state transition of alpha-HMM for multiple sequences of discrete continuous mixed series] The state transition probabilities are the same as in Equations 143 to 145. [Branch probability of alpha-HMM for discrete continuous mixed series of multiple arrays] The branching probability, that is, the grouping probability is expressed by the following equation as a modification of Equations 167 to 169.
  • the present embodiment also differs only in the sign and unknown parameter update formulas.
  • the hardware configuration shown in FIG. 8 is applied to the hardware configuration shown in FIG. 5 and the flowchart shown in FIG. 9 in the case of a plurality of arrangements, and detailed description thereof is omitted here.
  • the combination of sets of unknown parameters that are repeatedly updated is expressed by Formula 161
  • each update formula is expressed by Formula 130, Formula 128, and Formula 199 to 203 in the case of a single array as described above.
  • Equations 142 to 145 and Equations 204 to 210 are Equations 142 to 145 and Equations 204 to 210.
  • this invention is not limited to the said Example, It can change in the range which does not deviate from the meaning of this invention.
  • this invention is possible to increase a past value that goes back as in the modified example shown in the third embodiment, or to make a higher-order approximation.
  • the register 101 as a storage means for storing input observation data in time series, and the observation data are
  • the estimation unit 114 sets the value ⁇ of the HMM acceleration parameter.
  • Initial setting means 122 to perform, and the probability parameters of state transition probability a, branch probability c, mean value vector ⁇ , covariance matrix ⁇ , initial state probability ⁇ , and likelihood P (y
  • an update setting unit 124 for updating set an expected value N a and branched expected value N c of the state transition, the probability of immediately before the update set with update setting unit 124 and each Not only the waiting value but also each probability amount and each expected value shifted in time before that are used, and the observation data read from the register 101 and the speed-up parameter value ⁇ set by the initial setting means 122 are used.
  • each new probability amount and each expected value are calculated (refer to Formula 142 to Formula 145 and Formula 204 to Formula 210), and the calculation unit 126 converges the calculation. If the calculation is not converged, each new probability amount and each expected value calculated by the computing means 126 are set and updated by the update setting means, and if the computation is converged, the calculation means 126 is used. Determination means 128 for outputting each new probability amount as a final value.
  • the computing means 126 here calculates the value of the initial state probability ⁇ i
  • ⁇ l-1 ij is a subscript of a, l or l-1 is a subscript of ⁇
  • the above formulas 143 to 145 are used to calculate and
  • the update setting means 124 here sets the value of the acceleration parameter to 1 ⁇ ⁇ 3. By doing so, it is possible to reliably increase the speed of calculation of unknown parameters as compared to the conventional method by setting ⁇ to be larger than 1 in the range where ⁇ is less than 3 where repetitive calculations do not diverge.
  • the operational effects of the present embodiment described above are the initial setting step for executing the operation as the initial setting means 122, the update setting step for executing the operation as the update setting means 124, and the operation as the calculating means 126.
  • An HMM estimation method including a calculation step and a determination step for executing an operation as the determination unit 128 is also performed in exactly the same manner, and an HMM that causes a computer to function as such an estimation unit 114. Even if it is an estimation program of, it is demonstrated similarly.
  • the present embodiment can be applied to continuous series of observation data, it can be applied not only to speech recognition by a computer and robot motion recognition, but also to speech synthesis and robot motion generation.
  • the alpha-HMM algorithm applied in the present invention has very wide applicability such as speech recognition and synthesis, robot motion recognition and generation, neural information recognition, and life information sequence recognition.
  • the alpha-HMM algorithm proposed here is used to statistically handle temporal changes and stochastic fluctuations of the speech feature parameters as observation data. Can do.
  • parameters of the HMM are determined by inputting a plurality of speeches.
  • the probability can be calculated using an HMM whose parameters are determined from the speech to be recognized.
  • each probability amount is calculated based on discrete data.
  • speech synthesis or robot motion generation can be performed. It becomes possible.
  • the present invention recognizes genome sequences (for example, collects DNA information from one or more cancer patients, incorporates the data into the alpha-HMM algorithm proposed in the above embodiment, calculates the probability of developing cancer, etc.) ), Or robot motion recognition.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computational Linguistics (AREA)
  • Health & Medical Sciences (AREA)
  • Audiology, Speech & Language Pathology (AREA)
  • Human Computer Interaction (AREA)
  • Acoustics & Sound (AREA)
  • Multimedia (AREA)
  • Data Mining & Analysis (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Artificial Intelligence (AREA)
  • Compression, Expansion, Code Conversion, And Decoders (AREA)
  • Image Analysis (AREA)
  • Complex Calculations (AREA)

Abstract

　特殊な条件下以外であっても未知パラメータの計算が可能な隠れマルコフモデル推定アルゴリズムを得る。　推定手段１４は、隠れマルコフモデルの未知パラメータとして、状態遷移確率ａ，出力確率ｂ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび出力の期待値Ｎｂを更新設定する更新設定手段２４と、更新設定手段２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記レジスタ１から読み出した観測データと、初期設定手段２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する演算手段２６とを備える。

Description

隠れマルコフモデルの推定方法，推定装置および推定プログラム
　本発明は、あるデータ系列が与えられたときに、それを生成する隠れマルコフモデル（ＨＭＭ：Hidden Markov Model）を高速で推定するための隠れマルコフモデルの推定方法，推定装置および推定プログラムに関する。
　与えられたデータ系列から隠れマルコフモデルを推定するアルゴリズム（ＨＭＭ再推定アルゴリズム）は、1970年前後にほぼ同時に提案したエル．バウム（L.Baum）氏とエル．ウェルチ（L.Welch）氏の名前をとって、Baum-Welchアルゴリズムと呼ばれている。後に、これはマルコフ性の計算方式に特定したＥＭ（期待値最大化：Expectation-Maximization）アルゴリズムになっていることが分かっている。ＥＭアルゴリズムとは、観測データが現れる尤度（不完全データ尤度：incomplete-data likelihood）が最大となるような確率モデルを得るための手法であり、仮の完全データモデルに基づく対数尤度（log-likelihood）に対して、観測データに基づく条件付期待値を求めるステップ（E-step）と、その期待値の最大化を達成するようにモデルの更新を行なうステップ（M-step）とからなる二つの手順を、交互に繰り返すことで達成される。
　一方、ＥＭアルゴリズムは反復的であり、対数よりも有能な代理関数の使用が、汎用的で高速なアルゴリズム構造を提供すると期待される。そこで本願発明者の一人は、ＥＭアルゴリズムを部分集合として含むalpha-ＥＭアルゴリズムを、非特許文献１において提案している。これは、log（対数）関数の拡張版であるalpha-logを用いるものであり、ＨＭＭ再推定アルゴリズムについても、alpha-ＥＭアルゴリズム版（alpha-ＨＭＭ再推定アルゴリズム）があるものと予想されていた。
松山　泰男（Ｙ．Matsuyama），alpha-ＥＭアルゴリズム：alpha-log情報測定を利用した代理の尤度最大化法（The alpha-EM algorithm: Surrogate likelihood maximization using alpha-logarithmic information measures），IEEE情報理論会議（IEEE Trans. on Inform. Theory），第49巻（vol.49），pp.692-706，2003年
　しかし、上述したalpha-ＥＭアルゴリズムを用いた通常の導出法では、繰り返しの更新計算時に自分自身の値が必要となって計算できない、言い換えると更新計算する際に、その更新計算の結果と同時系列の項を使用して計算しなければならない、という自己撞着形となってしまい、特殊な条件下以外では、最終的に求めようとする確率モデルの計算ができない。そのため従来は、alpha-logを代理関数としたalpha-ＨＭＭ再推定アルゴリズムは不可能であるとみなされていた。
　そこで本発明は、上記問題点を解決して、特殊な条件下以外であっても未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることができる隠れマルコフモデルの推定方法，推定装置および推定プログラムを提供することを目的とする。
　本発明は、観測データを時系列に格納する記憶手段と、前記観測データがどのような確率モデルであるのかを、隠れマルコフモデルの未知パラメータを算出することで推定する推定手段とを備えた隠れマルコフモデルの推定装置において、前記推定手段は、前記隠れマルコフモデルの高速化パラメータを設定する初期設定手段と、前記隠れマルコフモデルの未知パラメータとして、状態遷移，出力，初期状態および尤度の各確率量と、状態遷移および出力の各期待値を更新設定する更新設定手段と、前記更新設定手段で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記記憶手段から読み出した観測データと、前記初期設定手段で設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算手段と、前記演算手段による計算の収束を判定し、計算が収束していなければ、前記演算手段で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、前記演算手段で計算した新たな各確率量を最終的な値として出力させる判定手段とを備えて構成される。
　この場合の演算手段は、前記新たな各確率量として、ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、後述する数６３の式で計算し、ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、数６４の式で計算し、状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、数６５の式で計算する構成とする。
　さらに更新設定手段は、前記高速化パラメータの値を１＜β＜３に設定する構成とする。
　本発明は、記憶手段に観測データを時系列に格納し、前記観測データがどのような確率モデルであるのかを、推定手段が隠れマルコフモデルの未知パラメータを算出することで推定する隠れマルコフモデルの推定方法において、前記未知パラメータの算出は、前記隠れマルコフモデルの高速化パラメータを設定する初期設定ステップと、前記隠れマルコフモデルの未知パラメータとして、状態遷移，出力，初期状態および尤度の各確率量と、状態遷移および出力の各期待値を更新設定する更新設定ステップと、前記更新設定ステップで更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記記憶手段から読み出した観測データと、前記初期設定ステップで設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算ステップと、前記演算ステップによる計算の収束を判定し、計算が収束していなければ、前記演算ステップで計算した新たな各確率量および各期待値を前記更新設定ステップで設定更新させ、計算が収束していれば、前記演算ステップで計算した新たな各確率量を最終的な値として出力させる判定ステップとにより行なわれる。
　この場合の演算ステップは、前記新たな各確率量として、ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を計算し、ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して計算し、状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して計算する。
　また、観測データが離散系列でない場合には、演算ステップは、アルゴリズムの繰り返し指標（インデックス）を示すｌ(エル）を付与した演算式により、状態遷移や分岐の期待値を利用して、さらに出力確率の代わりに分岐確率と、平均値ベクトルと、共分散行列とを計算する。
　また、観測データが複数本ある場合には、演算ステップは、その観測データが何本目の配列かを示すインデックスを付して、同様に計算する。
　さらに更新設定ステップは、前記高速化パラメータの値を１＜β＜３に設定する。
　本発明は、記憶手段に時系列に格納された観測データがどのような確率モデルであるのかを、隠れマルコフモデルの未知パラメータを算出することで推定する推定手段として、コンピュータを機能させる隠れマルコフモデルの推定プログラムにおいて、前記推定手段を、前記隠れマルコフモデルの高速化パラメータを設定する初期設定手段と、前記隠れマルコフモデルの未知パラメータとして、状態遷移，出力，初期状態および尤度の各確率量と、状態遷移および出力の各期待値を更新設定する更新設定手段と、前記更新設定手段で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記記憶手段から読み出した観測データと、前記初期設定手段で設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算手段と、前記演算手段による計算の収束を判定し、計算が収束していなければ、前記演算手段による新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、前記演算手段による新たな各確率量を最終的な値として確定させる判定手段として機能させている。
　この場合の演算手段は、前記新たな各確率量として、ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を計算し、ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して計算し、状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して計算する。
　また、観測データが離散系列でない場合には、演算手段は、アルゴリズムの繰り返し指標（インデックス）を示すｌを付与した演算式により、状態遷移や分岐の期待値を利用して、さらに出力確率の代わりに分岐確率と、平均値ベクトルと、共分散行列とを計算する。
　また、観測データが複数本ある場合には、演算手段は、その観測データが何本目の配列かを示すインデックスを付して、同様に計算する。
　さらに更新設定手段は、前記高速化パラメータの値を１＜β＜３に設定している。
　請求項１，４，１４の発明によれば、ＨＭＭの未知パラメータとして、状態遷移，出力，初期状態および尤度の各確率量と、状態遷移および出力の各期待値を算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることができる。
　請求項２，５～１２，１５～２２の発明によれば、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
　請求項３，１３，２３の発明によれば、繰り返しの計算が発散しない範囲で、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
本発明の第１実施例を示す推定装置の処理手順を示すフローチャートである。 同上、推定装置の構成を示すブロック図である。 同上、alpha-ＨＭＭアルゴリズムの収束速度について、繰り返し回数と尤度との相関関係を示すグラフである。 本発明の第２実施例を示す推定装置の処理手順を示すフローチャートである。 同上、推定装置の構成を示すブロック図である。 離散型アルファベットの場合のモデル図である。 連続型アルファベットの場合のモデル図である。 本発明の第３実施例、第５実施例、及び第７実施例の単一配列のデータの場合に対応する推定装置の処理手順を示すフローチャートである。 本発明の第４実施例、第６実施例、及び第７実施例の複数配列のデータの場合に対応する推定装置の処理手順を示すフローチャートである。
　１，１０１　レジスタ（記憶手段）
　１４，１１４　推定手段
　２２，１２２　初期設定手段
　２４，１２４　更新設定手段
　２６，１２６　演算手段
　２８，１２８　判定手段
　以下、添付図面を参照しながら、本発明における隠れマルコフモデルの推定方法，推定装置および推定プログラムの好ましい各実施例について説明する。
　先ず、本発明の具体的な実施例を提示する前に、当該実施例における独自のalpha-ＨＭＭ再推定アルゴリズムを導入するまでの理論的な経緯について、以下の数式を参照しながら説明する。
　一般に、隠れマルコフモデルによる学習アルゴリズムは、凸ダイバージェンスの最適化から始まって生成されている。まず、次の式に示すように、２つの確率密度p，qの間の凸ダイバージェンスについて考える。
Figure JPOXMLDOC01-appb-M000141
　ここで、ＹはＫ次元ユークリッド空間である。関数ｆ(r）は、(０，∞）において凸である。関数f（r）の双対関数g（r）は、次式を満たす必要がある。
Figure JPOXMLDOC01-appb-M000142
　上記数１における「≧」という記号は、ｐ＝ｑを殆どどこでも保持する場合に限り保持する。ｆ（1）の正規化が任意であるため、ここでは、Ｄｆ（ｐ・ｑ）とＤｇ（ｑ・ｐ）がｐとｑとの間の擬似距離として作用するように、ｆ（1）を次のように選定する。
Figure JPOXMLDOC01-appb-M000143
　留意すべきなのは、次の２つの式が、カルバック・ライブラー・ダイバージェンス（Kullback-Leibler divergence）を発生させる、ということである。もしｐとｑが確率量の関数であれば、数１における積分は総和に置き換えられる。これは、以後説明する積分の全てにおいて同様である。
Figure JPOXMLDOC01-appb-M000144
Figure JPOXMLDOC01-appb-M000145
　上記数１～数５の関係は、対数より有能であるものの、対数と類似した関数であることを意味している。連続して２回微分可能な関数のクラスについて考えると、数２は以下の等式を与える。
Figure JPOXMLDOC01-appb-M000146
　また、ｒ＝１付近では以下の式がそれぞれ成り立つ。
Figure JPOXMLDOC01-appb-M000147
Figure JPOXMLDOC01-appb-M000148
　ここで、o（1）は高次無限小の項を意味する。もし、乗算的に分離可能な凸関数を使用するのであれば、o（1）は不要であることに注目されたい。
Figure JPOXMLDOC01-appb-M000149
　以下の関数は、その一例である。
Figure JPOXMLDOC01-appb-M000150
　もし、上記の関数を、ｆ（1）＝０およびｆ”（1）＝１であるようにシフトして正規化すれば、ｋ＝－ｃ（１－ｃ）となる。前記数７のコアとなる項は以下の式を有している。
Figure JPOXMLDOC01-appb-M000151
　上式において、次の式は興味深い関数であり、これはパラメータｃによって凸性が調整できる単調関数である。
Figure JPOXMLDOC01-appb-M000152
これは、以下の式の関係が保たれることを指し示すのに重要である。
Figure JPOXMLDOC01-appb-M000153
　したがって、Ｌ(c)（r）はパラメータｃを有する法則化された対数と見なされる。つまり、ｃ-対数と呼ぶことができる。しかしながら、パラメータｃの増加はその数値を減少させる。そこで、ここでは以下の関係によって、法則化された対数をパラメータ化する。
Figure JPOXMLDOC01-appb-M000154
　このように、α-対数（alpha-log）は以下の式に示すような凸ダイバージェンスＤｆ（ｐ・ｑ）から始まることによって得られる。
Figure JPOXMLDOC01-appb-M000155
　ここでの対数は次式のようになる。
Figure JPOXMLDOC01-appb-M000156
　注目すべきは、次の式に示す凸関数が、ｒ=ｐ/ｑとして上記数１に適用されるとαダイバージェンスが発生し、これはalpha-ＥＭアルゴリズムにおいて重要な役割を果たす。
Figure JPOXMLDOC01-appb-M000157
　αダイバージェンスは、次の式であらわせる。
Figure JPOXMLDOC01-appb-M000158
　上式において、ｐ＝ｑが殆どどこでも適用される場合に限り、等式のゼロが達成される。
　次に、alpha-ＨＭＭ再推定アルゴリズムの原点となるalpha-ＥＭアルゴリズムの期待値と最大化の処理ステップについて説明する。ここではＰＹ｜ψ（ｙ｜ψ）を、ψによってパラメータ化されたＹにおける観測データｙの確率密度または確率量とする。集合Ｙは、不完全データの集合と見なされる。ｘ∈Ｘ（太字）は、未知または消失したデータを含む理想的観測の完全データあるいは増補データとする。そして、不完全データのpdf（probability density function：確率密度関数）あるいはpmf（probability mass function：確率量関数）は次式に示される。
Figure JPOXMLDOC01-appb-M000159
　ここでの積分範囲は、次式のようになり、pdfの積分はpmfの総和となる。
Figure JPOXMLDOC01-appb-M000160
　そして、条件付きのpdfまたはpmfは次式の通りとなる。
Figure JPOXMLDOC01-appb-M000161
　pdfまたはpmfの不完全なデータのalpha-log尤度比は、次式の通りである。
Figure JPOXMLDOC01-appb-M000162
　ここで、次の式におけるαダイバージェンスの計算は、alpha-ＥＭアルゴリズムの基本的な関係を示している。
Figure JPOXMLDOC01-appb-M000163
　その式は、以下の通りである。
Figure JPOXMLDOC01-appb-M000164
　したがって、最大化のための目的関数は以下の式のような定量化（quantity）によって表される。
Figure JPOXMLDOC01-appb-M000165
　上記説明によって、alpha-ＥＭアルゴリズムがalpha-ＨＭＭ再推定アルゴリズムの開始点と位置付けられることがわかる。
　alpha-ＥＭアルゴリズムでは、次のステップで処理が行なわれる。
　初期設定：最初のサイクル用のパラメータφを設定する。
　Ｅ-ステップ：φが与えられたときのψの関数として、前記数２５を計算する。
　Ｍ-ステップ：次の式にしたがって、更新パラメータを計算する。
Figure JPOXMLDOC01-appb-M000166
　Ｕ-ステップ：φをψ*に置換して、収束するまではＥ-ステップに戻る。 
もし、閉じた部分で正確な最大化が得られなければ、数２６に示す最大化が改善した計算に置き換わる。
　また、別なalpha-ＧＥＭアルゴリズムでは、次のステップで処理が行なわれる。
　初期設定：最初のサイクル用のパラメータφを設定する。
　Ｅ-ステップ：φが与えられたときのψの関数として、前記数２５を計算する。
　Ｍ-ステップ：次の式にしたがって、更新パラメータを計算する。
Figure JPOXMLDOC01-appb-M000167
　Ｕ-ステップ：φをψ＋に置換して、収束するまではＥ-ステップに戻る。alpha-ＨＭＭの近似されたバージョンが、ＧＥＭアルゴリズムと見なすことができる点に注目されたい。
　隠れマルコフモデル（ＨＭＭ）の公式化の問題は、一連のランダムな観測変数データを発生させる最良のマルコフモデルを推定することにある。観測データの変数列は、次の式であらわせる。この式におけるtはマルコフモデルを生成するためのデータの順序を表すものであり、繰り返しのインデックスとは意味が異なる。
Figure JPOXMLDOC01-appb-M000168
　ここでの「最良」とは、隠れマルコフモデルの尤度が最大化されることを意味する。隠れマルコフモデルを公式化するために、次に定義する各確率を推定する必要がある。
　（ａ）状態遷移確率
　状態がｉからｊに遷移する確率のことであり、次式のように定義される。
Figure JPOXMLDOC01-appb-M000169
　上記式において、ｓi∈Σは有限な範囲での総和Σの中の１つの状態を示すアルファベットの記号である。状態ｓkから状態ｓlへのつながりが無い場合、状態遷移確率ａkl＝０となる。これは、状態遷移における従前のトポロジーを反映している。
　（ｂ）出力確率
　状態ｓjにおいて出力ｙkが発生する確率のことであり、次式のように定義される。
Figure JPOXMLDOC01-appb-M000170
　ここで、上記式における各要素は以下の通りである。
Figure JPOXMLDOC01-appb-M000171
　（ｃ）初期状態確率
　最初の状態がｉである確率のことであり、次式のように定義される。
Figure JPOXMLDOC01-appb-M000172
　推定される上記（ａ）～（ｃ）の各要素は、以下のように示される。
Figure JPOXMLDOC01-appb-M000173
　つまり、ここで用いる記号θは、状態遷移確率，出力確率および初期状態確率を含ませたものである。
　一方、隠れマルコフモデル推定の問題に関し、前記数２８で示したランダムな観測データの変数列は、その値が次の式であらわされる。ここで、tはデータの順序を表すインデックスである。
Figure JPOXMLDOC01-appb-M000174
　与えられた変数列は、その変数列を発生できる最大尤度モデルθを推定する。この隠れマルコフモデル推定の問題は、以下の解釈をもたらす。
（ａ）数２８で示す観測データが、前述した不完全データＹとして与えられる。
（ｂ）次の式に示す状態遷移は、消失データＺである。
Figure JPOXMLDOC01-appb-M000175
（ｃ）完全データは、次式のように示される。
Figure JPOXMLDOC01-appb-M000176
　この確率は、次式で示される。
Figure JPOXMLDOC01-appb-M000177
（ｄ）前述の三重項θは、数３６に示す完全データの確率構造を提供する。
　この解釈は期待値最大化によって、実際にはalpha-ＥＭアルゴリズムによって、前記ＨＭＭ問題の公式化を可能にする。そして、以下の重要な定理が得られる。
　定理：ψとφを、２つのマルコフモデルの各パラメータセットとする。そして、α≦１に対して、次の数３８の式は、数３９という意味を持つ。
Figure JPOXMLDOC01-appb-M000178
Figure JPOXMLDOC01-appb-M000179
　ここで数３９の不等式は、数３８の不等式が成立する場合に成立するものである。
　証明：基本の方程式である数２４において、次式に示す不完全データのalpha-log尤度比を設定する。
Figure JPOXMLDOC01-appb-M000180
　数３５に示す状態Ｓ（太字）が消失データであるとすると、Ｘ＝（Ｓ，Ｙ）（各記号は何れも太字）であるので、前記数２５のＱ-関数は次式で算出できる。
Figure JPOXMLDOC01-appb-M000181
　前記数１８の特性と数２４の基本的な方程式によって、この定理が保たれる。
　数４１の特性は、数４０の形態で推定されるマルコフモデルψが、マルコフモデルφよりも高い確率で数列ｙ（太字）を生成することを意味する。すなわち、 数３３に示すマルコフモデルクラスでのalpha-ＥＭアルゴリズムの繰り返しは、必然的な局所最適性のコストで最良の隠れマルコフモデルを生成する。したがって、数３３に示すマルコフモデルのパラメータで数３８を繰り返すことが、ＨＭＭアルゴリズムの理論版となる。ここでは、α＝１の場合（後述する高速化パラメータの値が、β＝３である場合)は、収束限界であるために実用的ではない。
　alpha-ＥＭアルゴリズムとしての主な問題は、ソフトウェアとして実現できる具体的なコンピュータアルゴリズムをいかに提供するのかにある。数３８のalpha-ＨＭＭは理論上の形態であり、これはlog-ＥＭアルゴリズムと同様に一般的なＨＭＭである。ソフトウェアで実行可能なalpha-ＨＭＭとして、ここでは２つのバージョンを提示するが、第１のバージョンは正確であるものの、未来の情報の計算を必要とする。すなわち、得られたアルゴリズムはnon-causalである。一方、本発明で特に提案する第２のバージョンは、時間シフトと確率の近似を用いることによって、未来の情報を含めることなく計算が可能となる。なお、ここで提示するalpha-ＨＭＭアルゴリズムの生成において、全てはpmfであるので、積分が総和となる。以下、それぞれのバージョンについて説明する。
　（Ａ）第１の正確なバージョン
　・non-causalな更新
　正確なバージョンは、条件付き微分によって得られる。このバージョンで得られる更新がnon-causalであることは、既に述べている。この形態は次のバージョンのcausalな形態を得るための原点となるものである。
　次式に示すように、全ての更新式は、ｘ＝（ｓ，ｙ）（各記号は何れも太字）に対する前記Ｑ-関数の最大化から得られる。
Figure JPOXMLDOC01-appb-M000182
　ここで、ｔはalpha-ＨＭＭに対する繰り返しの指標（インデックス）である。そして、パラメータ集合は次式に示すようになる。なお、明細書中、小文字のπ、a、b等の記号は集合として表記する場合は大文字で示すようにしている。
Figure JPOXMLDOC01-appb-M000183
　・Πθt+1（Πは太字、t+1はθの添え字）の推定
　初期状態確率πｉ｜θt+1（πは太字、t+1はθの添え字）の更新は、次の修正されたＱ-関数における微分の最大化によって得られる。
Figure JPOXMLDOC01-appb-M000184
　この計算は、次式の微分を用いることによって達成できる。
Figure JPOXMLDOC01-appb-M000185
　そして、次式を用いることで、数４４のλを除去する。
Figure JPOXMLDOC01-appb-M000186
　すると、次の更新式を与える。
Figure JPOXMLDOC01-appb-M000187
　上記数４７は観測データｙ（太字）に対して計算できるように見えるが、左辺に含まれるθt+1が右辺に含まれている。これは、non-causalityの存在による自己撞着性を示している。しかしながら、これは次のバージョンの近似によって解決されることとなる。したがって、確率のための２つの更なる更新式を得ることができた。
　・Ａθt+1（Ａは太字、t+1はθの添え字）の推定
　状態遷移確率ａij｜θt+1（t+1はθの添え字）の更新は、同様に微分によって得られる。
Figure JPOXMLDOC01-appb-M000188
　ここで次式を用いることで、λを除去する。
Figure JPOXMLDOC01-appb-M000189
　すると、次の更新式を与える。
Figure JPOXMLDOC01-appb-M000190
　ここでのＮij（ｓ）（ｓは太字）は、状態ｓにおいてｉからｊに移行する状態遷移の数である。数５０も右辺にθt+1を含んでいる。これも次のバージョンで解決されることとなる。
　・Ｂθt+1（Ｂは太字、t+1はθの添え字）の推定
　出力確率ｂjk｜θt+1（t+1はθの添え字）の更新は、再度の微分によって得られる。
Figure JPOXMLDOC01-appb-M000191
　ここで次式を用いることで、λを除去する。
Figure JPOXMLDOC01-appb-M000192
　すると、次の更新式を与える。
Figure JPOXMLDOC01-appb-M000193
　ここでのＮｂjk（ｓ）（ｓは太字、jkはｂの添え字）は、数列ｓ（太字）によって生じる出力確率ｂjkに対する事象の数である。数５３も右辺にθt+1を含んでいる。これは次のバージョンで解決されることとなる。
　（Ｂ）第２のCausalな近似バージョン
　前述したように、第１のバージョンに示す一連の更新式（数４７，数５０，数５３）は、自己撞着性を解決するために近似を必要とする。
　・時間のシフト
　alpha-ＥＭアルゴリズムの収束は、適切な収束判定基準によって、次の式のようになることを意味する。
Figure JPOXMLDOC01-appb-M000194
　これは、以下の近似式を与える。ここで、o(1)は高位の無限小を意味する。
Figure JPOXMLDOC01-appb-M000195
　これは、繰り返し指標tをシフトすることによって、数４７，数５０，数５３に示す更新式の自己撞着性を解決することになる。
　・αのシフト
　前記時間のシフトは、尤度比の期待値のための確率環境が反転されることを意味する。したがって、元のパラメータαはαcausalに変換される。αとαcausalの関係は、以下の誘導によって得られる。ここで、数４７，数５０，数５３における尤度比の期待値の核となる部分を考慮すると、数４７の場合は以下の関係が成り立つ。
Figure JPOXMLDOC01-appb-M000196
　この誘導は尤度比の時系列を維持しながら、期待値が現在の環境Ｐ（ｓ｜ｙ，θt）（ｓ，ｙは太字）から、未来の環境Ｐ（ｓ｜ｙ，θt+1）（ｓ，ｙは太字）に変化することを意味する。上記数５６を見るとθｔ＋１の初期確率πiを左辺として計算するのに、第一式の右辺では同じ時系列のθｔ＋１に依存する関数（項）をしなければならないという問題点があり、このままでは計算できないと言う自己撞着性を含んでいた。そこで、まず第二式では第一式の分母分子を｛｝内の項のα乗でそれぞれ除算する。さらに第三式ではテーラー展開によって微小近似すると共に時間シフトの概念を導入しており、これによって自己撞着矛盾を解決している。言い換えると、右辺には左辺と同時系列のθｔ＋１に依存する関数（項）を用いない形になっている。この詳細な考え方を以下に示す。
　ａ）未来に関する尤度比を現在の環境で計算すると、次式のようになる。
Figure JPOXMLDOC01-appb-M000197
　ｂ）未来に関する尤度比を未来の環境で計算すると、次式のようになる。
Figure JPOXMLDOC01-appb-M000198
　したがって、時間シフト前後の関係は以下の通りである。
Figure JPOXMLDOC01-appb-M000199
　そこから得られる関係式は、次式のようになる。
Figure JPOXMLDOC01-appb-M000200
　したがって、α＝－１のlog-ＥＭとlog-ＨＭＭは、ここでのαcausal＝１の場合に相当する。なお、符号の使用を簡単化するために、以後、高速化パラメータの値となるαcausalを、次式のようにβとしてあらわす。すると、数６０の右側の条件式はβ≦３を示すことになり、従来技術であるlog-ＥＭとlog-ＨＭＭはβ＝１の場合に対応することがわかる。（ただし、β＝３だと後述の実験結果より発散してしまうことがわかっている。）
Figure JPOXMLDOC01-appb-M000201
　・テーラー展開
　前記全ての更新式（数４７，数５０，数５３）は、尤度比の計算に能力を必要とする。これは計算に時間がかかるものであり、Ｔが増加するに従って扱いにくくなる。テーラー展開は、alpha-logアルゴリズムにおける尤度比の長所を失うことなく、計算を簡単化することができるもので、尤度比は次式のように近似化される。
Figure JPOXMLDOC01-appb-M000202
　数５５と数６２を数４７，数５０，数５３に適用すると、過去の情報を利用するcausalで、かつ計算上効率的なalpha-ＨＭＭアルゴリズムが得られる。注目すべきは、これが数２７のalpha-ＧＥＭアルゴリズムに対応するということである。
　（Ｃ）ソフトウェアで実行可能なalpha-ＨＭＭアルゴリズム
　・Πθt+1（Πは太字、t+1はθの添え字）の推定
　高次の項ｏ（1）を放棄することによって数５５と数６２の近似式を結合し、数４７に適用すると以下の更新式が与えられる。
Figure JPOXMLDOC01-appb-M000203
　・Ａθt+1（Ａは太字、t+1はθの添え字）の推定
　数５５と数６２の近似式を適用すると、数５０は以下の更新式を与える。
Figure JPOXMLDOC01-appb-M000204
　ここで、Ｎａij｜θt（ijはａの添え字、tはθの添え字）は、条件付き確率Ｐ（ｓ｜ｙ，θt）（ｓ，ｙは太字）の下でのＮij（ｓ）（ｓは太字）の期待値であり、不完全データｙ（ｙは太字）によって生じる状態遷移のカウント数により推定される。
　・Ｂθt+1（Ｂは太字、t+1はθの添え字）の推定
　数５５と数６２の近似式を適用すると、数５３は以下の更新式を与える。
Figure JPOXMLDOC01-appb-M000205
　ここで、ＮｂjK｜θt（jkはｂの添え字、tはθの添え字）は、条件付き確率Ｐ（ｓ｜ｙ，θt）（ｓ，ｙは太字）の下でのＮｂjK（ｓ）（ｓは太字）の期待値であり、不完全データｙ（ｙは太字）によって生じる状態遷移のカウント数によって推定される。
　・従来のlog-ＨＭＭとの比較
　上記の数６３，数６４，数６５は、ソフトウェアとして実行可能である。これらの更新式は、以下に示す共通の特性を有する。
（１）β＝１すなわちα＝－１の場合は、log-ＥＭアルゴリズムから得られた従来のＨＭＭアルゴリズムと一致する。
（２）数６０と数６１との間の関係により、αが行なうのと同様に、パラメータβは更新用の現在と過去の情報の各要素を調整する。
（３）βに依存する、すなわちαに依存するすべての過去の更新値は、格納データを参照することによってのみ得ることができる。
（４）現在の項は、log-ＨＭＭアルゴリズムに対する前方－後方アルゴリズムを用いることで効率的に計算できる。
　上記特性（１）～（４）は、alpha-ＨＭＭアルゴリズムの各更新サイクルがlog-ＨＭＭアルゴリズムのオーバーヘッドより僅かに増えてだけであることを示唆している。したがって、設計パラメータα、あるいは高速化パラメータβが適切に選択されれば、本発明のalpha-ＨＭＭアルゴリズムは、従来のlog-ＨＭＭアルゴリズムよりも高速に収束することが予測される。この違いが生じる理由は、数６３、数６４、数６５の式において、従来のlog-ＥＭから得られたＨＭＭアルゴリズムに相当するβ＝１を代入すると分母分子の第二項が消去され、分母分子の第一項が約分されて尤度比が残らなかったのに対し、β＞1で設定すると分母分子の第二項が消去されず、単純に分母分子を約分できないために尤度比の影響を残すことができるようになったことだと考えられる。
　以下、上記理論に基づく好ましい実施例について説明する。
　実施例１、２では、観測データが離散的な場合において、観測データが単一配列（一本の観測データ）を有する場合、複数配列（Ｍ本の観測データ）を有する場合について、それぞれ説明する。
　実施例３、４では、観測データが連続的な場合において、観測データが単一配列（一本の観測データ）を有する場合、複数配列（Ｍ本の観測データ）を有する場合について、それぞれ説明する。
　実施例５、６では、観測データが半連続的な場合において、観測データが単一配列（一本の観測データ）を有する場合、複数配列（Ｍ本の観測データ）を有する場合について、それぞれ説明する。
　実施例７では、観測データが離散的かつ連続的な場合、すなわち離散系列と連続系列が混在している場合において、単一配列を有する場合と複数配列を有する場合をまとめて説明する。
　図１は、本発明の第１実施例において、上述したalpha-ＨＭＭ再推定アルゴリズムを実行可能にするプログラムの処理手順をあらわしたものである。
　同図において、１は合計でＴ個の観測データを格納する記憶手段としてのレジスタで、各観測データは時刻τが１から順に時系列に並んで格納される。ここでｙτは個々のデータ値を示し、ｙはＴ個のデータ値の集合を示しており、本実施例では配列が１本の観測データをレジスタ１に格納している。推定装置１０は、レジスタ１に記憶される一列のデータ値が、どのような確率構造（モデル）を有しているのかを、ＨＭＭの未知パラメータを算出する以下のステップＳ１～ステップＳ１０の各手順に従って推定解析するものである。
　ステップＳ１は、前述した高速化パラメータの値βを設定する部分である。従来のＨＭＭ推定アルゴリズムは、β＝１の場合に相当する。つまり、β＝１という特殊な条件下であれば、従来のalpha-ＨＭＭ再推定アルゴリズムであっても、上記確率モデルの計算は可能である。しかし、本実施例で提案する新たなalpha-ＨＭＭ再推定アルゴリズムでは、確率モデルの計算が可能な高速化パラメータの値βを１≦β＜３の範囲に拡張することができる。これは、従来のβ＝１を特例として含むものである。なお、βが大きな値であるほど、推定装置１０としての処理の高速性は増すが、収束性を保持するにはβ＜３でなければならない。
　ステップＳ２は、初期確率と収束判定値を決める部分である。これは、後述するステップＳ３～ステップＳ８の手順を繰り返す前に行なわれる。推定装置１０が最終的に算出しようとするＨＭＭの確率構造は、次式のようにあらわせる。これは、前述の数３７と等しい。
Figure JPOXMLDOC01-appb-M000206
　上記式において、πS0（０はｓの添え字）は最初の状態ｓ０における確率であり、ａSτ-1Sτ（τ－１およびτはｓの添え字）は時刻τ－１の状態ｓτ-1が時刻τの状態ｓτに移る確率であり、ｂSτ（ｙτ）（τはｓの添え字）は時刻τの状態ｓτに移行したときにｙτが出力される確率である。推定装置１０は、レジスタ１に格納された観測データを読み出して、その観測データが最も出現しやすくなる初期状態確率πと、状態遷移確率ａと、出力確率ｂとを推定するが、このステップＳ２では、次式に示すように、それらの確率π，ａ，ｂの初期値を決定する。
Figure JPOXMLDOC01-appb-M000207
　上式において、θ０は数４３にあるように、０回目の繰り返しにおける初期の確率π，ａ，ｂの組合せを示し、推定装置１０は、その条件で状態ｉから状態ｊに移る状態遷移確率ａij｜θ０（０はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂjk｜θ0（０はθの添え字）の値と、最初に状態ｉとなる初期状態確率πｉ｜θ0（０はθの添え字）の値とをそれぞれ決定する。
　またステップＳ２では、対数尤度に基づく収束範囲を決めることで、ステップＳ３～ステップＳ９の繰り返しを終了させるための収束判定値を決定する。この収束判定値はステップＳ８で用いられ、具体的には後述する数８３で示される。
　これらの値をステップＳ１，Ｓ２で決定すると、ステップＳ３の手順に移行して、確率量とカウント値の設定が行なわれる。ステップＳ３は、最初にステップＳ２で求めた初期確率を利用し、それ以降はステップＳ９で更新された確率量とカウント値を利用して、数６８に示す各確率値と、数６９に示す各カウント値を、実際に推定装置１０のメモリ（図示せず）に設定する部分である。
Figure JPOXMLDOC01-appb-M000208
Figure JPOXMLDOC01-appb-M000209
　次のステップＳ４は、レジスタ１から観測データを読み出して、ステップＳ３で設定された上式の確率量から、数７０と数７１に示す前向き確率を計算する部分である。なお、ここに示すαは確率値であり、前述したパラメータの値αとは異なる。ここでは、数７０に示す確率値αが計算され、そこから数７１に示す確率値（尤度）Ｐ（ｙ｜θt）が計算される。
Figure JPOXMLDOC01-appb-M000210
Figure JPOXMLDOC01-appb-M000211
　次のステップＳ５は、レジスタ１から観測データを読み出して、ステップＳ３で設定された上式の確率量から、数７２に示す前向き確率を計算する部分である。なお、ここに示すβは確率値であり、前述したパラメータの値βとは異なる。ここでは、数７２に示す確率値βが計算される。
Figure JPOXMLDOC01-appb-M000212
　上記前向き確率と後向き確率は、計算回数を減らすために、既存のＨＭＭ推定アルゴリズムにも組み込まれていたものである。
　続くステップＳ６は、ステップＳ４で計算された確率値αと、ステップＳ５で計算された確率値βを用い、レジスタ１から観測データを読み出して、数７３に示す状態遷移のカウント値と、数７４に示す出力のカウント値をそれぞれ計算する部分である。
Figure JPOXMLDOC01-appb-M000213
Figure JPOXMLDOC01-appb-M000214
　次のステップＳ７は、t+1回目に繰り返される確率π，ａ，ｂの組合せを条件として、状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値と、状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値とを、ステップＳ４で求めたt回目の確率Ｐ（ｙ｜θt）と、それよりも１回前の確率Ｐ（ｙ｜θt-1）を利用し、且つステップＳ５で求めたt回目のカウント値Ｎａij｜θt，Ｎｂjk｜θt（ijはａの添え字、jkはｂの添え字、tはθの添え字）と、それよりも１回前のカウント値Ｎａij｜θt-1，Ｎｂjk｜θt-1（ijはａの添え字、jkはｂの添え字、t-1はθの添え字）を利用して計算する部分であり、これは前記数６３，数６４，６５に対応している。
　このステップＳ７では、前記ステップＳ１で設定したパラメータの値βが用いられ、本アルゴリズムの主要な特徴部分となっている。特に、添え字ｔ+1の部分の確率を計算する上で、時間シフトした添え字tや添え字t-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したことが、その効果として挙げられる。
　本実施例では、配列が１本の観測データをレジスタ１に格納しているが、その場合のステップＳ７における確率用の更新式は、次式のようにあらわせる。
Figure JPOXMLDOC01-appb-M000215
Figure JPOXMLDOC01-appb-M000216
Figure JPOXMLDOC01-appb-M000217
　ステップＳ８は、計算の収束を判定する部分である。ここでは、直前のステップＳ７で計算された新たな確率量に基づく尤度Ｐを用い、その尤度Ｐと前記ステップＳ２における収束判定値との比較により、ステップＳ３～ステップＳ９の計算が収束したか否かを判定する。近似的には、繰り返しの回数を指定して、指定回数に達したら、計算が収束したと判定してもよい。
　ステップＳ８において計算が収束していないと判定された場合、ステップＳ９に移行して手順が繰り返される。ステップＳ９は、計算された確率量と事象のカウント値を、繰り返し計算のために更新する部分である。本実施例では、次式に示すように、２つの過去値をシフトさせる。
Figure JPOXMLDOC01-appb-M000218
　ここでは、繰り返しの回数が１つ増えることにより、t回目の確率量およびカウント値がt-1回目の確率量およびカウント値に更新され、t+1回目の確率量およびカウント値がt回目の確率量およびカウント値に更新される。なお、ステップＳ３とステップＳ９を一纏めにして、更新された確率量およびカウント値をそのまま次の繰り返しのために、推定装置１０のメモリに設定記憶させてもよい。
　一方、ステップＳ８において計算が収束していると判定された場合、ステップＳ１０に移行して、次式に示す計算された三組の確率量を用いたＨＭＭを採用する。
Figure JPOXMLDOC01-appb-M000219
　上記数７９で計算した各値を用いて、推定装置１０は数６６に示すＨＭＭの確率構造を出力することができる。
　図２は、図１に示すアルゴリズムを実現する推定装置１０の構成を示している。同図において、推定装置１０は前記レジスタ１を内蔵する入力手段１２と、レジスタ１に格納した観測データを適宜読み出し、上記ステップＳ１～ステップＳ１０の手順を実行して数６６に示すＨＭＭの確率構造を推定する推定手段１４と、この推定手段１４で得たＨＭＭの確率構造を出力する出力手段１６とにより概ね構成される。本実施例の入力手段１２は、１本の配列を有する観測データをレジスタ１に格納できる構造となっているが、複数本の配列を有する観測データを取り扱う推定装置１００については、次の実施例で詳しく説明する。
　推定装置１０のハードウェア構成は、例えば演算処理部としてＣＰＵを備えたコンピュータで実現することができる。その場合、図１に示す推定装置１０の処理手順を実行するプログラムが、メモリなどの記録媒体に記憶される。当該プログラムをどこに記憶するのかは限定せず、例えば通信手段を介してプログラムがコンピュータにダウンロードされる構成であってもよい。
　推定手段１４は、前記ステップＳ１における高速化パラメータの値βや、ステップＳ２における初期の確率値および収束判定の条件を設定する初期設定手段２２と、ステップＳ３における確率量およびカウント値の設定や、ステップＳ９における確率量およびカウント値の更新を行なう更新設定手段２４と、ステップＳ４における前向き確率の計算や、ステップＳ５における後ろ向き確率の計算を行なうと共に、そこからステップＳ６における状態遷移のカウント値と出力のカウント値をそれぞれ計算し、さらにステップＳ７における新たな確率量の計算を行なう演算手段２６と、ステップＳ８における計算の収束判定を行ない、計算が収束していなければ、前記更新設定手段２４による確率量およびカウント値の更新を行なわせる一方で、計算が収束していれば、ステップＳ１０において、演算手段２６が直前に計算した新たな確率量を、ＨＭＭの最終的なパラメータ値として確定させる判定手段２８とを備えている。
　また推定装置１０には、前記高速化パラメータの値や、初期の確率値および収束判定の条件や、更新設定される確率量およびカウント値の値を読み書き可能に格納するメモリ３０の他に、必要に応じて高速化パラメータの値や収束判定の条件を操作入力するキーボードやマウスなどの操作手段３２が、推定手段１４に接続して設けられる。これらの各装置構成により、上述したステップＳ１～ステップＳ１０の手順が実行される。
　次に、本実施例における試験用の観測データでの実験結果を説明する。
　（Ａ）生成されたデータを用いた速度の評価
　観測データの集合が強い凸性でなければ（そして強い凹面でもなければ）、全ての最適化アルゴリズムは局所的に最適となる。存在するデータの殆ど全ては、そのような性質を持っている。従来のあるいはlog-ＨＭＭアルゴリズムは例外的ではないし、alpha-ＨＭＭアルゴリズムも例外的ではない。したがって、そうしたＨＭＭアルゴリズムを適用する前に、発生のメカニズムが分かっている人工的なデータ集合を生成することが必要になってくる。
　ここでは、数３３に示す数列をレジスタ１に格納する入力データとし、そのパラメータを次式のように指定する。
Figure JPOXMLDOC01-appb-M000220
　このマルコフ連鎖を使用し、[0, 1]上の一様な乱数によって、1000個のサンプルを有するデータ集合を生成した。完全な推定には、以下の困難があることを先に述べる。
（ａ）各データ集合の殆どが、局所的最適を与える対象である。
（ｂ）たとえ局所的最適が避けられたとしても、正確な推定には無限の長さの数列が必要である。
（ｃ）コンピュータが生成した乱数は単なる擬似的な乱数にすぎない。
　したがって、以下が満たされれば満足するものとする。
（１）数３３に示すパラメータ集合の数値が、数８０からずれないこと
（２）alpha-ＨＭＭアルゴリズムの収束したパラメータ集合がlog-ＨＭＭのそれと近いこと
　ステップＳ２における初期値は、次のように設定する。これにより、初期状態を除いて事前の情報は必要としないものとする。
Figure JPOXMLDOC01-appb-M000221
　以下の理由のため初期状態は固定した。
（ａ）分布された初期確率は、不安定な局所最適性を得やすい。これは、alpha-ＨＭＭアルゴリズムの各収束速度を比較するのに不要な曖昧さをもたらす。
（ｂ）ＨＭＭ推定の後の認識問題において、一つの初期状態はビタビ（Viterbi）アルゴリズムによって選択される。
　図３は、対数尤度（log-likelihood）、すなわち繰り返し回数を考慮した尤度Ｐ（ｙ｜θt）に関する収束の傾向を示している。ここでは、数６０および数６１を満たす多様な値で実験した。
　この図から以下のことがわかる。
（１）β=αcausalが増加するにつれて対数尤度の立ち上がりが早くなる。すなわち、収束するまでの繰り返し回数が少なくなる。
（２）alpha-ＥＭアルゴリズムの限界(α＝１)は、alpha-ＨＭＭアルゴリズムではβ=αcausal＝３に相当するが、この場合、発散してしまう。
　さらに２つの視点で、alpha-ＨＭＭアルゴリズムの性能を比較する必要がある。一つは局所最適の比較であり、もう一つはＣＰＵ時間である。局所最適の比較では、次式に示すように、θに対する全パラメータと共に、それに付随した対数尤度の値を一覧にする必要がある。
Figure JPOXMLDOC01-appb-M000222
　上記の実験結果の一覧を見るとβを1～2.75の範囲で変化させても対数尤度（ＬＬ）の値がほとんど変化せずに正確さを保っていることがわかる。この実験ではβ≧３に設定した実験は発散する結果になるために省略しているが、３に近い2.75に設定しても実験結果にそれほど影響なく高速化できることがわかる。ただし、このような実験結果の収束限界となるβの値は入力データによって若干異なる。
　次は、収束するまでのＣＰＵ時間の比較である。この比較はソフトウェアの実行手段とハードウェアの手法に依存する。ここでは、alpha-ＨＭＭアルゴリズムを標準的なPOSIX環境におけるＣ‐コードとして実行する。これらのコードは標準的なシングルコアプロセッサで走らせている。前記ステップＳ８で行われる収束は、以下の判定基準で測定した。
Figure JPOXMLDOC01-appb-M000223
　上式において、Ｐnewは今回の計算で得た確率Ｐの値であり、Ｐoldは前回の計算で得た確率Ｐの値である。この収束判定基準は更新前後で確率の値がほとんど変化していなければ収束したと判断するものである。この収束判定式結果は次の表に示す通りであり、ここでは繰り返し回数とＣＰＵ時間による収束の比較を示している。
Figure JPOXMLDOC01-appb-T000224
　局所最適性とＣＰＵ時間とを合せて比較すると、以下のことがわかる。
（ａ）収束する各対数尤度は非常に接近しているものの、何れも局所最大値が存在する。
（ｂ）対数のまたは従来のＨＭＭアルゴリズムは、常に最良の対数尤度をもたらす訳ではない。log-ＨＭＭアルゴリズムよりも良好な局所最大値が、alpha-ＨＭＭアルゴリズムで得られた。
（ｃ）改良したalpha-ＨＭＭアルゴリズムによって、繰り返し回数とＣＰＵ時間の高速化が成し遂げられたのは明らかである。
（ｃ）繰り返し回数の改善が、直接ＣＰＵ時間に影響している。これは前記数６３，数６４および数６５に示すように、予め計算されてメモリに格納される過去の情報を利用しているからであり、従来のlog-ＨＭＭアルゴリズムに比べて、ソフトウェアの複雑さを増やしてはいない。
（ｄ）一方、αcausal＝β＝３すなわちα＝１の場合は、発散する。これは、alpha-ＥＭアルゴリズムの許容範囲と一致する。したがって、数６２による近似は許容される。αcausal＝β＝１すなわちα＝－１の場合において、数６３，数６４および数６５はlog-ＨＭＭアルゴリズムの正確な形態であることが重要である。この意味において、数６３，数６４および数６５の更新式はlog-ＨＭＭアルゴリズムの拡張版ということができる。
　前記図３を参照すると、従来のalpha-ＨＭＭアルゴリズムで唯一ステップＳ６からステップＳ７への計算が可能であったβ＝１の結果に対して、今回提案するalpha-ＨＭＭアルゴリズムで計算が可能になったβ＝2.75の結果では、少ない繰り返し回数で推定装置１０の処理が収束していることが判る。具体的には、上記表１において、β＝１の場合は繰り返し回数が263回で、ＣＰＵの使用時間が0.253秒となっているが、β＝2.75の場合は繰り返し回数が70回、またＣＰＵの使用時間が0.068秒に減少する。表１の最右列にあるのは、β＝１のときの繰り返し回数とＣＰＵの使用時間を１としたときの速度比を計算したもので、例えばβ＝2.75に設定すれば繰り返し回数は3.76倍、またＣＰＵの使用時間は3.72倍に高速化する。
　以上の実験結果から、本実施例では従来よりも高速なＨＭＭアルゴリズムが示された。本推定方法は、基になっているalpha -ＥＭアルゴリズムを反映しており、改良したalpha-ＨＭＭアルゴリズムと呼ぶことができる。改良したalpha-ＨＭＭアルゴリズムは、従来のalpha-ＨＭＭアルゴリズムまたはlog-ＨＭＭアルゴリズムより優れている。計算上の複雑さの増加が非常に少ないので、繰り返し回数の抑制が直接的に期待したＣＰＵの高速化を実現した。
　また、改良したalpha-ＨＭＭアルゴリズムで採用する数６３，数６４および数６５の更新式は、メモリに格納した過去の情報を追加的に必要とするだけである。したがって、既存のコードからのソフトウェアのバージョンアップは難しくない。この場合、alpha-対数の曲率を制御するαcausal＝α＋２が設計パラメータとなる。
　さらに、より少ないサンプルで未来の情報を利用する方法について、確認を行なった。この方法は、テーラー展開と共にnon-causalな形態を使用したことにある。更新式の一例は、次式のように表せる。
Figure JPOXMLDOC01-appb-M000225
　ここで、θ＾ｔ＋１とθ＾ｔ（＾はθの上部に記される）は、より少ないサンプル（例えば1000個のサンプルから200個を取り出す）から統計量が推定されることを示している。すなわち、 サンプルの一部が未来を考えるのに使われる。この方法は数６３，数６４および数６５の更新式よりも大きなオーバーヘッドを必要とするが、高速化の利益の方がこの余計な重荷よりも勝っている。もし、入力源がよく混合されている、あるいはエルゴードであると演繹的知識をユーザーが持っていれば、この方法を利用できる。実験では、高速化が数６３，数６４および数６５の方法よりも悪くなかった。
　以上のように本実施例では、入力される観測データを時系列に格納する記憶手段としてのレジスタ１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１４とを備えたＨＭＭの推定装置１０において、推定手段１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，出力確率ｂ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび出力の期待値Ｎｂを更新設定する更新設定手段２４と、更新設定手段２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数６８，数６９を参照）と共に、前記レジスタ１から読み出した観測データと、初期設定手段２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数６３，数６４，数６５，数７５，数７６，数７７を参照）演算手段２６と、演算手段２６による計算の収束を判定し、計算が収束していなければ、演算手段２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段２６で計算した新たな各確率量を最終的な値として出力させる判定手段２８とを備えている。
　このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，出力確率ｂ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび出力の期待値Ｎｂを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
　また、ここでの演算手段２６は、ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、前記数６３の式で計算し、ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、前記数６４の式で計算し、状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、前記数６５の式で計算する構成となっている。
　そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
　さらに、ここでの更新設定手段２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
　なお上述した本実施例の作用効果は、初期設定手段２２としての動作を実行する初期設定ステップと、更新設定手段２４としての動作を実行する更新設定ステップと、演算手段２６としての動作を実行する演算ステップと、判定手段２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。
　図４は、本発明の第２実施例におけるプログラムの処理手順をあらわしたものである。
　同図において、１０１は合計でＴ個の観測データを格納する記憶装置としてのレジスタで、各観測データのそれぞれは、時刻τが１から順に時系列に並んで格納される。ここでｙτは個々のデータ値を示し、ｙはＴ個のデータ値の集合を示し、ｎは何本目の配列かをあらわすインデックスを示しており、本実施例では複数の配列を有するＭ（＝２以上の整数）本の観測データをレジスタ１０１に格納している。推定装置１００は、レジスタ１０１に記憶されるＭ列のデータ値が、どのような確率構造（モデル）を有しているのかを、以下のステップＳ１１～ステップＳ２０の各手順に従って推定解析するものである。
　ステップＳ１１は、前述した高速化パラメータの値βを設定する部分である。これは上記第１実施例のステップＳ１に相当するもので、高速化パラメータの値βは、従来のβ＝１を特例として含む１≦β＜３の範囲に拡張される。なお、βが大きな値であるほど、推定装置１００としての処理の高速性は増すが、収束性を保持するにはβ＜３でなければならない。
　ステップＳ１２は、初期確率と収束判定値を決める部分である。これは、前記ステップＳ２に相当するもので、後述するステップＳ１３～ステップＳ１８の手順を繰り返す前に行なわれる。推定装置１００は、レジスタ１０１に格納された観測データを読み出して、その観測データが最も出現しやすくなる初期状態確率πと、状態遷移確率ａと、出力確率ｂとを推定するが、このステップＳ１２では、次式に示すように、それらの確率π，ａ，ｂの初期値を決定する。
Figure JPOXMLDOC01-appb-M000226
　上式において、θ０は数４３にあるように、０回目の繰り返しにおける初期の確率π，ａ，ｂの組合せを示し、推定装置１００は、その条件で状態ｉから状態ｊに移る状態遷移確率ａij｜θ０（０はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂjk｜θ0（０はθの添え字）の値と、最初に状態ｉとなる初期状態確率πｉ｜θ0（０はθの添え字）の値とをそれぞれ決定する。
　またステップＳ１２では、対数尤度に基づく収束範囲を決めることで、ステップＳ１３～ステップＳ１９の繰り返しを終了させるための収束判定値を決定する。この収束判定値はステップＳ１８で用いられ、具体的には前記数８３で示される。
　これらの値をステップＳ１１，Ｓ１２で決定すると、ステップＳ１３の手順に移行して、確率量とカウント値の設定が行なわれる。ステップＳ１３は、最初にステップＳ１２で求めた初期確率を利用し、それ以降はステップＳ１９で更新された確率量とカウント値を利用して、数８６に示す各確率値と、数８７に示す各カウント値を、実際に推定装置１００のメモリ（図示せず）に設定する部分である。
Figure JPOXMLDOC01-appb-M000227
Figure JPOXMLDOC01-appb-M000228
　次のステップＳ１４は、レジスタ１０１から観測データを読み出して、ステップＳ３で設定された上式の確率量から、数８８と数８９に示す前向き確率を計算する部分である。なお、ここに示すαは確率値であり、前述したパラメータの値αとは異なる。ここでは、数８８に示す確率値αが計算され、そこから数８９に示す確率値（尤度）Ｐ（ｙ(n)｜θt）が計算される。
Figure JPOXMLDOC01-appb-M000229
Figure JPOXMLDOC01-appb-M000230
　次のステップＳ１５は、レジスタ１０１から観測データを読み出して、ステップＳ１３で設定された上式の確率量から、数９０に示す前向き確率を計算する部分である。なお、ここに示すβは確率値であり、前述したパラメータの値βとは異なる。ここでは、数９０に示す確率値βが計算される。
Figure JPOXMLDOC01-appb-M000231
　上記前向き確率と後向き確率は、計算回数を減らすために、既存のＨＭＭ推定アルゴリズムにも組み込まれていたものである。
　続くステップＳ１６は、ステップＳ１４で計算された確率値αと、ステップＳ５で計算された確率値βを用い、レジスタ１０１から観測データを読み出して、数９１に示す状態遷移のカウント値と、数９２に示す出力のカウント値をそれぞれ計算する部分である。
Figure JPOXMLDOC01-appb-M000232
Figure JPOXMLDOC01-appb-M000233
　次のステップＳ１７は、t+1回目に繰り返される確率π，ａ，ｂの組合せを条件として、状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値と、状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値とを、ステップＳ１４で求めたt回目の確率Ｐ（ｙ(n)｜θt）と、それよりも１回前の確率Ｐ（ｙ(n)｜θt-1）を利用し、且つステップＳ５で求めたt回目のカウント値Ｎ(n) ａjk｜θt，Ｎ(n) ｂjk｜θt（jkはａの添え字、tはθの添え字）と、それよりも１回前のカウント値Ｎ(n) ａjk｜θt-1，Ｎ(n) ｂjk｜θt-1（jkはａの添え字、t-1はθの添え字）を利用して計算する部分であり、これは前記数６３，数６４，６５に対応している。
　このステップＳ１７では、前記ステップＳ１１で設定したパラメータの値βが用いられ、本アルゴリズムの主要な特徴部分となっている。特に、添え字ｔ+1の部分の確率を計算する上で、時間シフトした添え字tや添え字t-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したことが、その効果として挙げられる。
　本実施例では、配列がＭ本の観測データをレジスタ１０１に格納しているが、その場合のステップＳ１７における確率用の更新式は、次式のようにあらわせる。
Figure JPOXMLDOC01-appb-M000234
Figure JPOXMLDOC01-appb-M000235
Figure JPOXMLDOC01-appb-M000236
　ステップＳ１８は、計算の収束を判定する部分である。ここでは、直前のステップＳ１７で計算された新たな確率量に基づく尤度Ｐを用い、その尤度Ｐと前記ステップＳ１２における収束判定値との比較により、ステップＳ１３～ステップＳ１９の計算が収束したか否かを判定する。近似的には、繰り返しの回数を指定して、指定回数に達したら、計算が収束したと判定してもよい。
　ステップＳ１８において計算が収束していないと判定された場合、ステップＳ１９に移行して手順が繰り返される。ステップＳ１９は、計算された確率量と事象のカウント値を、繰り返し計算のために更新する部分である。本実施例では、次式に示すように、２つの過去値をシフトさせる。
Figure JPOXMLDOC01-appb-M000237
　ここでは、繰り返しの回数が１つ増えることにより、t回目の確率量およびカウント値がt-1回目の確率量およびカウント値に更新され、t+1回目の確率量およびカウント値がt回目の確率量およびカウント値に更新される。
　一方、ステップＳ１８において計算が収束していると判定された場合、ステップＳ２０に移行して、次式に示す計算された三組の確率量を用いたＨＭＭを採用する。
Figure JPOXMLDOC01-appb-M000238
　上記数９７で計算した各値を用いて、推定装置１０は数６６に示すＨＭＭの確率構造を出力することができる。
　図５は、図４に示すアルゴリズムを実現する推定装置１００の構成を示している。同図において、推定装置１００は前記レジスタ１０１を内蔵する入力手段１１２と、レジスタ１０１に格納した観測データを適宜読み出し、上記ステップＳ１１～ステップＳ２０の手順を実行して数６６に示すＨＭＭの確率構造を推定する推定手段１１４と、この推定手段１１４で得たＨＭＭの確率構造を出力する出力手段１１６とにより概ね構成される。本実施例の入力手段１１２は、Ｍ本の配列を有する観測データをレジスタ１０１に格納できる構造となっている。推定装置１００のハードウェア構成については、前記第１実施例の推定装置１０と同様であるため、ここでは説明を省略する。
　推定手段１１４は、初期設定手段１２２と、更新設定手段１２４と、演算手段１２６と、判定手段１２８とを備えている。これらは前記第１実施例の初期設定手段２２，更新設定手段２４，演算手段２６および判定手段２８にそれぞれ対応するもので、取り扱う観測データがＭ本の配列になった以外は、第１実施例と同様に機能する。
　また推定装置１１０には、前記高速化パラメータの値や、初期の確率値および収束判定の条件や、更新設定される確率量およびカウント値の値を読み書き可能に格納するメモリ１３０の他に、必要に応じて高速化パラメータの値や収束判定の条件を操作入力するキーボードやマウスなどの操作手段１３２が、推定手段１１４に接続して設けられる。これらの各装置構成により、上述したステップＳ１１～ステップＳ２０の手順が実行される。
　以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１０１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１１４とを備えたＨＭＭの推定装置１００において、推定手段１１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段１２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，出力確率ｂ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび出力の期待値Ｎｂを更新設定する更新設定手段１２４と、更新設定手段１２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数６８，数６９を参照）と共に、前記レジスタ１０１から読み出した観測データと、初期設定手段１２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数６３，数６４，数６５，数９３，数９４，数９５を参照）演算手段１２６と、演算手段１２６による計算の収束を判定し、計算が収束していなければ、演算手段１２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段１２６で計算した新たな各確率量を最終的な値として出力させる判定手段１２８とを備えている。
　このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，出力確率ｂ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび出力の期待値Ｎｂを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
　また、ここでの演算手段１２６は、ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、前記数６３の式で計算し、ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、前記数６４の式で計算し、状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、前記数６５の式で計算する構成となっている。
　そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
　さらに、ここでの更新設定手段１２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
　なお上述した本実施例の作用効果は、初期設定手段１２２としての動作を実行する初期設定ステップと、更新設定手段１２４としての動作を実行する更新設定ステップと、演算手段１２６としての動作を実行する演算ステップと、判定手段１２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。
　上記の実施例１、２では、観測データが離散的な場合について述べた。本実施例では、観測データが単一配列で連続的な場合の実施形態について説明する。
　この観測データが連続的な場合とは、観測データのデータ列の各々の値が波のように変動しているような場合である。
　具体的には、図６の観測データが離散的な場合のモデル図における出力確率ｂik，ｂjkが、図７の観測データが連続的な場合のモデル図のように確率密度関数である分岐確率ｃik，ｃjkに置き換わるというものである。
　（Ａ）前記第１実施例と前記第２実施例の式の変形
　ここで、本実施例の具体的な内容に入る前に、離散的な場合について上記実施例と同様な説明を重複して行なう。
　これは、実施例１、２と実施例３～８との間において、数式の簡略化や符号の変更等が入っており、その点を踏まえて整合性の取れた説明をするという便宜上の理由のためである。なお、上記実施例１、２と実質的な差異はない。
　まず、ソースデータとマルコフモデルの各種パラメータについて説明すると、HMMによるモデル化において、次式のようなソースデータ列が与えられる。ここで、tはデータの順序を表すインデックスである。
Figure JPOXMLDOC01-appb-M000239
　各ｙtは、スカラーあるいはベクトルである。上式は、単一配列の場合であるが、M本の複数配列で与えられるならば、ソースデータ列は次式の通りである。記号｛｝は集合を示している。
Figure JPOXMLDOC01-appb-M000240
　まず、単一配列では、数９８の観測データｙ（ｙは太字）が与えられる。この場合のHMMの課題は、次式の最尤推定法（Maximum Likelihood Estimation：MLE）の認識において、最良なモデルを見つけることである。
Figure JPOXMLDOC01-appb-M000241
　ここで、S（Sは太字）は、次式のような状態遷移系列を意味している。
Figure JPOXMLDOC01-appb-M000242
　ランダムな変数とそれらの値のために、小文字のｙ（ｙは太字）およびｓ（ｓは太字）は同様に用いられる（通常、ランダム変数は大文字で示される）。
　各確率は以下の数１０２～数１０５の通りである。
　・初期状態確率
Figure JPOXMLDOC01-appb-M000243
　・状態遷移確率
Figure JPOXMLDOC01-appb-M000244
　ここで、もし、状態ｓｉから状態ｓｊへのつながりが無ければ、状態ｉから状態ｊに移る状態遷移確率はａij＝０となる。これは、状態遷移のための前のトポロジーを反映している。
　・出力確率
Figure JPOXMLDOC01-appb-M000245
　ここで、上式の各要素は次式の通りである。
Figure JPOXMLDOC01-appb-M000246
　上記数１０２～数１０５に示した確率の集合をまとめてθとして次式に示す。
Figure JPOXMLDOC01-appb-M000247
　そして、確率的なデータ構造について以下のように解釈している。不完全なデータはｙ（ｙは太字）であり、推定される消失データはｓ（ｓは太字）であり、完全なデータは、数１００の確率量と数１０６のパラメータを有しており、次式に示されるｘ（ｘは太字）として定義する。
Figure JPOXMLDOC01-appb-M000248
・Alpha-EMアルゴリズム
　不完全なデータ、消失データ、完全なデータによるHMMの解釈は、EMアルゴリズムに匹敵する。本紙の目的は、新たなHMM推定アルゴリズムを見つけることであるため、alpha-EMアルゴリズムから始まる道筋を示している。
　観測データｙ（ｙは太字）は推定すべきパラメータの全体を表すψによってパラメータ化されるため、Ｐy|ψ（ｙ｜ψ）（ｙは太字）を確率密度あるいは確率量とする。ｘ∈Ｘ（ｘとＸは太字）を、消失データを含む理想的な観測結果である完全なデータ、あるいは拡張されたデータとする。そして、不完全なデータ確率密度関数(pdf)、あるいは確率量関数(pmf)は次式で示される。
Figure JPOXMLDOC01-appb-M000249
　ここで、積分する範囲を次式とする。
Figure JPOXMLDOC01-appb-M000250
　pdfの積分は、pmfの総和となる。そして、条件付きのpdfあるいはpmfは次式である。
Figure JPOXMLDOC01-appb-M000251
　alpha-EMアルゴリズムでは、次式のようなalpha対数が用いられる。
Figure JPOXMLDOC01-appb-M000252
　ここで、α＝－１の場合は対数、すなわち、次式である。
Figure JPOXMLDOC01-appb-M000253
　alpha-EMアルゴリズムは、alpha対数に関して、不完全なデータの尤度比を考慮する必要があるため、次式のようになる。
Figure JPOXMLDOC01-appb-M000254
　ここで、φとψは、繰返し最大化ステップにおける数１０６に対して、古いモデルのパラメータと新しいモデルのパラメータを示している。それから、alpha-EMアルゴリズムの基本的な方程式が次式として得られる。
Figure JPOXMLDOC01-appb-M000255
　ここで、D(α)は、常に負とならない２つの条件付きの確率Pｘ|ｙ,φ（x|y,φ）（ｘとｙは太字）とPｘ|ｙ,ψ（x|y,ψ）（ｘとｙは太字）の間のalphaダイバージェンスである。上式において、次式のＱ関数が重要である。Ｅは定量化（quantity）によって表される最大化のための目的関数を示している。
Figure JPOXMLDOC01-appb-M000256
　数１１４のために、もし、このＱ関数が正ならば、a＜1の範囲で、数１１４の左辺である不完全なデータのalpha対数尤度比も正となる。したがって、alpha-EMアルゴリズムとその変化形であるalpha-GEMアルゴリズムについて以下に示す。
　[Alpha-EMアルゴリズム]
初期化: 数１０６の初期値を選択し、φとして用いる。
E-Step: 数１１５の計算を実行する。
M-step: 更新パラメータを次式で算出する。
　なお、"arg max"とは最大値を与える変数を意味する記号である。すなわち、次式は、Ｑ（α）（ｘ|ｙ、φ）（ψ｜φ）（ｘとｙは太字であり、（ｘ|ｙ、φ）はＱの下付き）を変数ψの関数と考えたとき、その最大値を与える変数の値をψ*とする、という意味である。
Figure JPOXMLDOC01-appb-M000257
U-step: φをψ*によって置き換え、収束することを確認する。収束しない場合は、E-stepに戻って繰り返しが繰返される。
　[Alpha-GEMアルゴリズム]
　これは、上記M-stepを次式のようにＱ関数を非負とするψ+の算出に置き換えたアルゴリズムである。
Figure JPOXMLDOC01-appb-M000258
　まず、alpha-HMMアルゴリズムの近似バージョンがalpha-GEMアルゴリズムであることに着目すべきである。数１１４を有するalpha-EMアルゴリズムの長所によるalpha-HMM推定アルゴリズムの基本的な特性について説明する。
　完全なデータをｘ＝（s，ｙ）（ｘ、ｓ、ｙは太字）とし、ｓ（ｓは太字）を消失データとし、ｙ（ｙは太字）を不完全なデータとし、数１１５は次式に相当するものとする。
Figure JPOXMLDOC01-appb-M000259
　もし、上式の値が負にならなければ、以下の不等式が成り立つ。
Figure JPOXMLDOC01-appb-M000260
　ここで、alpha-HMMの抽象バージョンに留意すべきである。
　・単一配列のalpha-HMM 
Non-Causalな更新方程式
　離散型のアルファベット系列ｙ（ｙは太字）の場合において、最大化されるQ関数は次式で示される。なお、ｌはアルゴリズムの繰返しの指標（インデックス）である。一方、上記文章中に出てくるｔはデータの順番に対応する指標（インデックス）である。したがって、ｌとｔは区別されるものである。
Figure JPOXMLDOC01-appb-M000261
　まず、状態遷移確率aijの更新方程式について説明する。状態遷移確率aij は、その更新後に確率量とする必要があるので、ラグランジュの未定係数法（Lagrange multiplier）を使用する必要がある。したがって、次式のように最大値を算出するための微分を行なう。
Figure JPOXMLDOC01-appb-M000262
　そうすると、上式から次式が得られる。
Figure JPOXMLDOC01-appb-M000263
　ここで、Ｎij(ｓ)（ｓは太字）は、iからjへの状態遷移の数と位置である。次のステップでは、ソフトウェアで実行可能なアルゴリズムを目指すため、以下の問題を解決する必要がある。
(a)Non-causalityの存在: 右辺にθl+1(θは太字）が含まれている。このままでは、自己撞着性のため計算できない。
(b)右辺の算出には、状態数がＮでデータ数がＴのとき、Ｏ(ＮT)すなわちＮのＴ乗のオーダーの演算が必要である。
　上記の問題は、次のセクションで解決されることとなるが、その前に、出力確率bjk|θl+1（l+1はθの添え字）と初期状態確率πi|θl+1（l+1はθの添え字）の２つの更なる更新式について次式に示す。なお、Ｎjk(ｓ)（ｓは太字）は、状態ｓt＝ｊでの出力yt＝kの発生回数である。
Figure JPOXMLDOC01-appb-M000264
Figure JPOXMLDOC01-appb-M000265
Causal近似と拡張系列: 離散出力の場合
　更新方程式である数１２２～数１２４のコアとなる部分は、繰返し指標（インデックス）のシフトによるCausal近似によって、次式のように変換することができる。
Figure JPOXMLDOC01-appb-M000266
　したがって、P（ｙ|θl）＝P（ｙ|θl-1）＋o(1)(yは太字）の領域で次の等式を得る。
Figure JPOXMLDOC01-appb-M000267
　ここで、αcausalという表記は、以後、βと表記する。
　いま、数１２５が算出可能であるが、演算の複雑度のために、もう一つの近似が必要となる。このため、我々は次式の系列拡張を用いている。
Figure JPOXMLDOC01-appb-M000268
　そして、causal近似式である数１２５の適用と数１２７の系列拡張は、以下に示す遷移確率の更新方程式を与える。なお、これらの確率の更新式である上記数１２２～数１２４は、それぞれ上記実施例１の数７５～７７に対応するものである。
　単一配列の離散型alpha-HMMの遷移確率aij|θl+1(l+1は、θの添え字)は、次式である。
Figure JPOXMLDOC01-appb-M000269
　ここで、以下の各特性について理解しておくことが重要である。
[特性１]
　β＝αcausal＝１の場合は、従来のlog-HMMの方法に帰趨する。
[特性２]
　数１２８の分子は、現在と過去の更新項の重み付けられた総和である。分母も同様である。
[特性３]
　数１２８の第２行及び第３行は、これらの確率計算に伴う複雑度を抑えた従来の前方-後方アルゴリズムの方法に匹敵する。
[特性４]
　alpha-HMMに唯一、追加する必要があるものは、θl-1での更新項を記憶することである。これは、CPU時間に直接現れているように繰返しを減らすことを意味する。実際の実験結果でも、この予想通りになっている。
　次に、残りの２つの更新方程式を示す。
　単一配列の離散型alpha-HMMの出力確率ｂjk|θl+1(l+1は、θの添え字)は、次式である。
Figure JPOXMLDOC01-appb-M000270
　単一配列の離散型alpha-HMMの初期状態確率πi|θl+1(l+1は、θの添え字)は、次式である。
Figure JPOXMLDOC01-appb-M000271
（変形例）
　以下、近似式の変形例について説明する。なお、以下の変形例は本発明の他の実施例でも同様に適用可能なものである。上述した隠れマルコフモデルの推定法では、有限個の過去情報を用いて各確率量をすべて次のような計算に基づいて計算している。
Figure JPOXMLDOC01-appb-M000272
　この高速化パラメータβは、利用者が選ぶデザインパラメータであり、現在値による繰り返し値に対し、1回前の繰り返し値を重みづけする効果を有している。この過去の値に重みづけを行って加算することは、有限の過去にまでさかのぼって行うこともできる。すなわち、ετ≧０をデザインパラメータとして設定し、次式のように置換してもよい。
Figure JPOXMLDOC01-appb-M000273
Figure JPOXMLDOC01-appb-M000274
そうすると、数１３１は次式となる。
Figure JPOXMLDOC01-appb-M000275
　次に展開による近似式を用いる場合について説明すると、上記数１３１や数１３４は、次式のように展開できる（例えばテーラー展開）。
Figure JPOXMLDOC01-appb-M000276
　ここで、o(1)は高位の無限小を意味する記号であり、上式の右辺第一項は，従来の隠れマルコフモデル推定アルゴリズムに相当している。そして、モーメンタム項に相当する右辺第二項は、展開により近似される高速化項の意味をもつ。
　なお、Ｔ＝１かつε１＝Ｐ（ｙ｜θl-1）／Ｐ（ｙ｜θl）（ｙはそれぞれ太字），すなわち重みを用いた場合は、重みを用いていない場合に比べて、収束速度はわずかに遅くなるものの、β＝３付近での収束性能の安定性が良くなる。換言すると、図３のような性能曲線において、計算が収束する立ち上がり部分では、高速化に伴い、波状の上下の変動が生じやすくなるが、その変動が少なくなる。
　・複数配列のalpha-HMM
　もし、HMMの予め設計されたトポロジーは、エルゴードなものであれば、単一の長いトレーニング系列ｙ（ｙは太字）は、十分である。もし、選択されたトポロジーが、吸収状態を有しているならば、複数のトレーニング系列を使うことが望まれる。このような複数系列に対するalpha-HMM推定の更新式は、上記の単一系列に対する方法を利用する形で得ることができる。
　離散的なシンボルの場合において、Ｓ（Ｓは太字）を次式のようにＭ個の状態遷移系列の集合とする。
Figure JPOXMLDOC01-appb-M000277
すると、複数配列のＱ関数は次式である。
Figure JPOXMLDOC01-appb-M000278
　ここで、Ｐは、次式に示されるMarkov過程の確率である。
Figure JPOXMLDOC01-appb-M000279
　上式において、初期状態確率πs0(n)(0はｓの下付きの添え字であり、（ｎ）は上付きの添え字である)、状態遷移確率ast-1(n)st(n) (t-1とtは、ｓの下付きの添え字であり、（ｎ）は上付きの添え字である)、及び出力確率ｂst(n)(yt (n)) (tは、ｓの下付きの添え字であり、（ｎ）は上付きの添え字である)の形態が、系列指標（インデックス）ｎから独立している点に着目すべきである。
　初期状態確率πi|θl+1(l+1は、θの添え字)によって、初期状態の更新方程式の誘導は、Ｑ関数の数１３７の微分から始まる。上述の単一配列との違いは、上式からも明らかなように微分がｎ回現れることである。そして、次式のようにnon-causalな方程式が得られる。
Figure JPOXMLDOC01-appb-M000280
　ここで、
Figure JPOXMLDOC01-appb-M000281
であり、上式の分母の関数ｆ(ｎ)は次式である。
Figure JPOXMLDOC01-appb-M000282
　そして、causalの繰返し指標（インデックス）のシフトと系列の拡張は、次の更新方程式を与える。
Figure JPOXMLDOC01-appb-M000283
　同様に、状態遷移と出力の更新方程式が得られる。
Figure JPOXMLDOC01-appb-M000284
　ここで、
Figure JPOXMLDOC01-appb-M000285
及び
Figure JPOXMLDOC01-appb-M000286
　また、次式の出力確率が得られる。
Figure JPOXMLDOC01-appb-M000287
　ここで、上式で省略された各項はそれぞれ次式となる。
Figure JPOXMLDOC01-appb-M000288
Figure JPOXMLDOC01-appb-M000289
　上記数１４７のΣの下にある記号「ｔ：ｙt＝ｋ」の意味は、ｙt＝ｋが成り立つｔのみを対象にする、という意味である。すなわち、上記数１４７の場合、ｙ＝｛ｙｔ｝T（ｔ＝１）（左辺のyは太字であり、ｔ＝１は｛ｙｔ｝の下付きの添え字）において、ｙt＝ｋが成り立つｔのみについて加算する、という意味である。
　（Ｂ）連続系列の単一配列
　ここから、連続的な出力系列の場合の本実施例の具体的な説明を行なう。
　もし、出力系列ｙ＝{ｙt}Ｔ（ｔ＝１）（出力系列ｙは太字であり、ｔ＝１は｛ｙt｝の下付き）が、連続的な多変数の観測結果として現れれば、数１２８～数１３０と似ているが、少しだけ異なる更新方程式が得られる。この場合、ｙtは太字の文字で示されていないが、適切な次元のユークリッド空間におけるベクトルである。
　このような連続的なアルファベットの場合、最尤推定法（MLE）の問題は、以下の尤度を最大化することである。
Figure JPOXMLDOC01-appb-M000290
　ここで、状態stで、出力確率ｂｓtｋt（ｙt）（ｓtｋtのtは、それぞれｓ，ｋの下付きの添え字）において、kt番目の枝へ遷移する確率を特定する分岐確率Ｃｓtｋt（ｓtｋtのtは、それぞれｓ，ｋの下付きの添え字）は、ｙtの確率密度関数である。我々は、これをガウス密度（Gaussian density）とみなしている。そして、状態ｊで状態ｋが出力される出力確率ｂｊｋ（ｙt）は次式となる。
Figure JPOXMLDOC01-appb-M000291
　ここで、μｊｋは平均値ベクトルであり、Σｊｋは共分散行列であり、総和と混同してはいけない。平均値ベクトルとは、分岐確率の確率密度関数の平均値を示している（以後、同様）。そして、状態jでの出力確率密度関数(pdf)は、ｂｊ（ｙt）がpdfとなるように次式となる。
Figure JPOXMLDOC01-appb-M000292
　ここで、上式の記号について説明する。
　まず、離散系列の場合、ｋを出力記号の種類を表すインデックスとすると、出力ykのとる値は有限種類であるため、ｂj(yk)はｂｊｋと略記できる。しかし、連続系列の場合は、離散値ではなく、連続値であるため、そのような表記はできない。そのため、時間を表すインデックスをｔとして、yｔはそのまま記載する必要がある。
　また、N（ｙｔ；μｊｋ，Σｊｋ）は、平均値ベクトルがμｊｋであり、Σｊｋを共分散行列とする多次元正規確率密度関数（多次元ガウス確率密度関数）で、ｙｔを分布の変数とするという意味である。
　実際、このようなガウス混合モデルは、「L. A Liporace, “Maximum likelihood estimation for multivariate observations of Markov sources” IEEE Trans. IT, vol. 28, pp. 729-734, 1982.」、及び「B.-H. Juang, “Maximum-likelihood estimation for mixture multivariate stochastic observations of Markov chains,” AT & T Tech. J., vol. 64, pp. 1235-1245, 1985.」等にも記載されているモデルであり、log-HMMでもここまでが導出可能な場合にとどまっている。
　図６は、離散的なアルファベットの場合を図示している。また、図７は、ガウス混合モデル（Gaussian mixture model）(bottom)の場合、すなわち連続的なアルファベットの場合を図示している。これらの図面を参照すると、図７中の分岐確率ｃｊｋの矢印は、図６中の出力確率ｂｊｋの矢印と対応していることがわかる。　数１４９と数１００を参照すると、数１４９の分岐確率ｃｓtｋt（ｓtｋtのtは、それぞれｓ，ｋの下付きの添え字）と出力確率ｂｓtｋt（ｙt）（ｓtｋtのtは、それぞれｓ，ｋの下付きの添え字）の積が、数１００の出力確率ｂｓt（ｙt）（ｓtのtは、ｓの下付きの添え字）に対応しているように見えるかもしれない。しかし、数１４９において離散値シンボルとみなされるのは、出力確率ｂｓtｋt（ｙt）（ｓtｋtのtは、それぞれｓ，ｋの下付きの添え字）を除いた部分であるため、図６、図７の対応関係との矛盾は生じない。
　混合確率の場合のため、消失データはｓ（ｓは太字）とｃ（ｃは太字）である。そのため、Ｑ関数は次式である。
Figure JPOXMLDOC01-appb-M000293
　数１４９、数１５２から明らかなように、初期確率と状態遷移確率の更新方程式は、それぞれ数１２８、数１３０と同様である。分岐確率ｃｊｋの更新方程式は、状態遷移確率aijの場合と同様にラグランジュの未定係数法（Lagrange multiplier）によって得られる。
Figure JPOXMLDOC01-appb-M000294
　そして、右辺は、状態遷移確率の場合と同様にcausalで計算可能となり、次式のように変形される。
Figure JPOXMLDOC01-appb-M000295
　次式は、平均値ベクトルμｊｋの更新方程式である。μｊｋについて数１５２の直接微分として、次のnon-causalの方程式が得られる。
Figure JPOXMLDOC01-appb-M000296
　そして、繰返し指標（インデックス）のシフト、系列の拡張、総和の変更により、次式が得られる。
Figure JPOXMLDOC01-appb-M000297
　数１５６が、θl-1の過去情報が十分に利用できる点を示していることに着目すべきである。
　共分散行列の更新には、行列微分が必要である。Ｑ関数の数１５２を共分散行列の逆行列Σ-1 ｊｋ（ｊｋはΣの添え字）について微分することにより、次のnon-causalな方程式が得られる。
Figure JPOXMLDOC01-appb-M000298
　そして、繰返し指標（インデックス）のシフト、系列の拡張、および総和の変更により、次の更新式が得られる。
Figure JPOXMLDOC01-appb-M000299
　ここで、上式の各項は以下の通りである。
Figure JPOXMLDOC01-appb-M000300
Figure JPOXMLDOC01-appb-M000301
　共分散行列の更新式である数１５８には、効果的に十分に利用できる過去情報の形態を有していることに着目すべきである。なお、記号Σ-1 ｊｋ｜θl;1とは、l+1回目の状態ｊで分岐ｋに移行したときの出力の共分散行列を表している（以後、同様）。
　ここで、ガウス混合（Gaussian mixture）alpha-HMM、すなわち単一系列の連続型alpha-HMMの更新方法について簡単に述べる。
[単一系列の連続型alpha-HMMの初期状態確率]
　更新式は、数１３０である。
[単一系列の連続型alpha-HMMの状態遷移確率]
　更新式は、数１２８である。
[単一系列の連続型alpha-HMMの分岐確率]
　更新式は、数１５４である。
[単一系列の連続型alpha-HMMの平均値ベクトル]
　更新式は、数１５６である。
[単一系列の連続型alpha-HMMの共分散行列]
　更新式は、数１５８であり、その要素は、数１５９、数１６０である。
　ここで、θl-1のすべての情報をメモリに記憶する点を改めて強調することが重要である。θlによって指標（インデックス）を付された項の算出は、log-HMMと同等である。
　図８は、本発明の第３実施例におけるプログラムの処理手順をあらわしたものである。このフローチャート全体の流れは、前記第１実施例の場合とほぼ同様である。ただし、本実施例では、観測データｙが、連続的な多変数の観測結果として現れる単一配列の連続系列データである点で異なる。また、初期値として設定され、繰返し計算される未知パラメータには、初期状態確率、状態遷移確率、分岐確率、平均値ベクトル、共分散行列があるため、各ステップにおける計算式も異なる。
　そのため、置き換わる数式の対応関係について説明する。まず、繰返し更新される未知パラメータの集合の組み合わせを次式とする。
Figure JPOXMLDOC01-appb-M000302
推定装置１０が最終的に算出しようとするHMMの確率構造は、数６６の代わりに数１４９となる。ステップＳ３２で決定する未知パラメータの初期値は、数６７の代わりに次式となる。
Figure JPOXMLDOC01-appb-M000303
　ステップＳ３２における収束判定値の決定方法は、数８３と同じである。
　ステップＳ３３では、最初にステップＳ３２で決定した未知パラメータの初期値を利用し、それ以降はステップＳ３９で更新された未知パラメータとカウント値を利用して、次の数１６３に示す各確率値と、数１６４に示す各カウント値を、実際に推定装置１０のメモリ（図示せず）に設定する。なお、前記実施例１、２と異なり、指標（インデックス）は、データの順番のｔではなく、アルゴリズムの繰り返しのｌを付して設定する。
Figure JPOXMLDOC01-appb-M000304
Figure JPOXMLDOC01-appb-M000305
　その後のステップＳ３４～Ｓ３６の演算動作は、前記実施例１の数７０～数７４の一部の符号を置き換えて適用することによって行う。具体的には、出力確率ｂｊ．ｙτ+1（τ＋１はｙの下付き添え字）を分岐確率ｃｊ．ｙτ+1（τ＋１はｙの下付添え字）と置き換え、指標（インデックス）であるｔをｌと置き換える。
　ステップＳ３７では、l+1回目に繰り返される確率π，a，ｃと、平均値ベクトルμ，及び共分散行列Σの組み合わせを条件として、状態iとなる初期状態確率
初期状態確率πｉ｜θl+1（l+1はθの添え字）の値と、状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂjk（ｙt）において、分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃｊｋ|l+1（l+1はθの添え字）の値と、平均値ベクトルμjｋ|θl+1（l+1はθの添え字）の値と，共分散行列Σjｋ|θl+1（l+1はθの添え字）の値と，ステップＳ３４で求めたl回目の確率Ｐ（ｙ｜θl）と、それよりも１回前の確率Ｐ（ｙ｜θl-1）を利用し、且つステップＳ３５で求めたl回目のカウント値Ｎａij｜θl，Ｎcjk｜θl（ijはａの添え字、jkはcの添え字、lはθの添え字）と、それよりも１回前のカウント値Ｎａij｜θl-1，Ｎcjk｜θl-1（ijはａの添え字、jkはｃの添え字、l-1はθの添え字）を利用して計算する部分であり、これは前記数１２８、数１３０、数１５４、数１５６、数１５８～数１６０の各更新式にそれぞれ対応している。
　このステップＳ３７においても高速化パラメータの値βが用いられ、添え字l+1の部分の確率を計算する上で、時間シフトした添え字lや添え字l-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したこと等は、前記実施例１と同様にその効果として挙げられる。
　本実施例では、配列が１本の観測データをレジスタ１に格納しているが、その場合のステップＳ３７における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数１３０、状態遷移確率の更新式は数１２８、分岐確率の更新式は数１５４、平均値ベクトルの更新式は数１５６、共分散行列の更新式は数１５８～数１６０である。
　ステップＳ３８では、前記実施例１と同様に、前記ステップ３７で計算された新たな確率量に基づく尤度Ｐを用いて、その尤度ＰとステップＳ３２における収束判定値との比較により、ステップＳ３３～Ｓ３９の計算が収束したか否かを判定する。収束してなければ、ステップＳ３９に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように２つの過去値をシフトさせる。
Figure JPOXMLDOC01-appb-M000306
　ここでは、繰り返しの回数が１つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。なお、ステップＳ３３とステップＳ３９を一纏めにして、更新された確率量およびカウント値をそのまま次の繰り返しのために、推定装置１０のメモリに設定記憶させてもよい。
　一方、ステップＳ３８において計算が収束していると判定された場合、ステップＳ４０に移行して、次式に示す計算された５組の確率量を用いたＨＭＭを採用する。
Figure JPOXMLDOC01-appb-M000307
　上記数１６６で計算した各値を用いて、推定装置１０は数１４９に示すＨＭＭの確率構造を出力することができる。
　図８に示すアルゴリズムを実現する推定装置の構成は、図２に示した前記第１実施例の推定装置１０と同様であるため、ここでは説明を省略する。
　以上のように本実施例では、入力される観測データを時系列に格納する記憶手段としてのレジスタ１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１４とを備えたＨＭＭの推定装置１０において、推定手段１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを更新設定する更新設定手段２４と、更新設定手段２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数１６３，数１６４を参照）と共に、前記レジスタ１から読み出した観測データと、初期設定手段２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１２８，数１３０，数１５４，数１５６，数１５８～数１６０を参照）演算手段２６と、演算手段２６による計算の収束を判定し、計算が収束していなければ、演算手段２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段２６で計算した新たな各確率量を最終的な値として出力させる判定手段２８とを備えている。
　このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
　また、ここでの演算手段２６は、l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、前記数１３０の式で計算し、l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、前記数１２８の式で計算し、状態ｊで状態ｋが出力される出力確率ｂjk｜θl+1（l+1はθの添え字）の際に分岐ｋに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1の値と、平均値ベクトルμjk｜θl+1（l+1はθの添え字）の値と、共分散行列Σjk｜θl+1の値とを、l回目およびl-1回目における前記出力の期待値Ｎcjk｜θlおよびＮcjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、前記数１５４、数１５６、数１５８～数１６０の式で計算する構成となっている。
　そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
　さらに、ここでの更新設定手段２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
　なお上述した本実施例の作用効果は、初期設定手段２２としての動作を実行する初期設定ステップと、更新設定手段２４としての動作を実行する更新設定ステップと、演算手段２６としての動作を実行する演算ステップと、判定手段２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。
　また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。
　本実施例では、連続的なシンボルで複数配列のデータ系列の場合の実施形態について説明する。
　本実施例においても、上記実施例と同様に、繰返し更新される未知パラメータの集合の組み合わせを数１６１とし、連続的なシンボルのalpha-HMMの更新方程式が得られる。初期状態確率と状態遷移確率の更新方程式は、それぞれ数１４２～数１４５と同様である。しかし、出力の更新方程式の集合は、離散的なシンボルの場合と異なる。我々は、分岐確率、平均値ベクトル、及び共分散行列の更新方程式を必要としている。
　分岐確率については、初期状態確率と状態遷移確率の場合と同様にラグランジュの未定係数法を用いて算出可能である。そして、次式の更新方程式が得られる。
Figure JPOXMLDOC01-appb-M000308
　ここで、各要素は次式である。
Figure JPOXMLDOC01-appb-M000309
Figure JPOXMLDOC01-appb-M000310
　Ｑ関数の数１３７での直接ベクトル微分であるμｊｋ｜θl+1（l+1はθの下付きの添え字）が適用される。その更新方程式は次式である。
Figure JPOXMLDOC01-appb-M000311
　ここで、
Figure JPOXMLDOC01-appb-M000312
　数１７０、数１７１と同様に、共分散行列の更新方程式は、Σ-1 ｊｋ｜θl（ｊｋ｜θlはΣの下付き添え字であり、さらにlはθの下付き添え字）について行列微分を用いることによって得られる。
Figure JPOXMLDOC01-appb-M000313
　ここで、Ｇ(n) θl（θlはＧの下付き添え字）は次式である。
Figure JPOXMLDOC01-appb-M000314
　図９は、本発明の第４実施例におけるプログラムの処理手順をあらわしたものである。このフローチャート全体の流れは、前記第３実施例の場合とほぼ同様である。ただし、本実施例では、観測データｙが、連続的な多変数の観測結果として現れる複数配列の連続系列データである点で異なる。
　置き換わる数式の対応関係について説明する。まず、推定装置１００が最終的に算出しようとするHMMの確率構造は、数６６の代わりに数１４９となる。ステップＳ４２で決定する未知パラメータの初期値は、数１６２と同じである。
　ステップＳ４２における収束判定値の決定方法は、数８３と同じである。
　ステップＳ４３では、最初にステップＳ４２で決定した未知パラメータの初期値を利用し、それ以降はステップＳ４９で更新された未知パラメータとカウント値を利用して、数１７４に示す各確率値と、数１７５に示す各カウント値を、実際に推定装置１００のメモリ（図示せず）に設定する。なお、指標（インデックス）は、データの順番のｔではなく、アルゴリズムの繰り返しのｌを付して設定する。
Figure JPOXMLDOC01-appb-M000315
Figure JPOXMLDOC01-appb-M000316
　その後のステップＳ４４～Ｓ４６の演算動作は、前記実施例２の数８８～数９２の一部の符号を置き換えて適用することによって行う。具体的には、出力確率ｂｊ．ｙτ+1（τ＋１はｙの下付添え字）を分岐確率ｃｊ．ｙτ+1（τ＋１はｙの下付添え字）と置き換え、指標（インデックス）であるｔをｌと置き換える。
　ステップＳ４７では、l+1回目に繰り返される確率π，a，ｃと、平均値ベクトルμ，及び共分散行列Σの組み合わせを条件として、状態iとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値と、状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂjk（ｙt）において、分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃｊｋ|l+1（l+1はθの添え字）の値と、平均値ベクトルμjｋ|θl+1（l+1はθの添え字）の値と，共分散行列Σjｋ|θl+1（l+1はθの添え字）の値とを、ステップＳ３４で求めたl回目の確率Ｐ（ｙ（n）｜θl）と、それよりも１回前の確率Ｐ（ｙ（n）｜θl-1）を利用し、且つステップＳ３５で求めたl回目のカウント値Ｎ（n） ａij｜θl，Ｎ（n） cjk｜θl（ijはａの添え字、jkはcの添え字、lはθの添え字）と、それよりも１回前のカウント値Ｎ（n） ａij｜θl-1，Ｎ（n） cjk｜θl-1（ijはａの添え字、jkはｃの添え字、l-1はθの添え字）を利用して計算する部分であり、これは前記数１４２～数１４５、数１６７～数１７３に対応している。
　このステップＳ４７においても高速化パラメータの値βが用いられ、添え字l+1の部分の確率を計算する上で、時間シフトした添え字lや添え字l-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したこと等は、前記実施例２と同様にその効果として挙げられる。
　本実施例では、配列がＭ本の観測データをレジスタ１０１に格納しているが、その場合のステップＳ４７における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数１４２、状態遷移確率の更新式は数１４３～数１４５、分岐確率の更新式は数１６７～数１６９、平均値ベクトルの更新式は数１７０、数１７１、共分散行列の更新式は数１７２、数１７３である。
　ステップＳ４８では、前記実施例２と同様に、前記ステップＳ４７で計算された新たな確率量に基づく尤度Ｐを用いて、その尤度ＰとステップＳ４２における収束判定値との比較により、ステップＳ４３～Ｓ４９の計算が収束したか否かを判定する。収束してなければ、ステップＳ４９に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように２つの過去値をシフトさせる。
Figure JPOXMLDOC01-appb-M000317
　ここでは、繰り返しの回数が１つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。
　一方、ステップＳ４８において計算が収束していると判定された場合、ステップＳ５０に移行して、次式に示す計算された５組の確率量を用いたＨＭＭを採用する。
Figure JPOXMLDOC01-appb-M000318
　上記数１７７で計算した各値を用いて、推定装置１０１は数１４９に示すＨＭＭの確率構造を出力することができる。
　図９に示すアルゴリズムを実現する推定装置の構成は、図５に示した前記第２実施例の推定装置１０１と同様であるため、ここでは説明を省略する。
　以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１０１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１１４とを備えたＨＭＭの推定装置１００において、推定手段１１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段１２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを更新設定する更新設定手段１２４と、更新設定手段１２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数１７４，数１７５を参照）と共に、前記レジスタ１０１から読み出した観測データと、初期設定手段１２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１４２～数１４５，数１６７～数１７３を参照）演算手段１２６と、演算手段１２６による計算の収束を判定し、計算が収束していなければ、演算手段１２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段１２６で計算した新たな各確率量を最終的な値として出力させる判定手段１２８とを備えている。
　このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
　また、ここでの演算手段１２６は、l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、前記数１４２の式で計算し、l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（t+1はθの添え字）の値を、l回目およびl-1回目における状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、前記数１４３～数１４５の式で計算し、状態ｊで状態ｋが出力される出力確率ｂjk｜θl+1（l+1はθの添え字）の際に分岐ｋの枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1の値と、平均値ベクトルμjk｜θl+1（l+1はθの添え字）の値と、共分散行列Σjk｜θl+1の値とを、l回目およびl-1回目における前記分岐の期待値Ｎcjk｜θlおよびＮcjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、前記数１６７～数１７３の式で計算する構成となっている。
　そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
　さらに、ここでの更新設定手段１２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
　なお上述した本実施例の作用効果は、初期設定手段１２２としての動作を実行する初期設定ステップと、更新設定手段１２４としての動作を実行する更新設定ステップと、演算手段１２６としての動作を実行する演算ステップと、判定手段１２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。
また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。
　本実施例では、半連続的なシンボルで単一配列のデータ系列の場合の実施形態について説明する。
　図７において、ガウス混合（Gaussian mixture）alpha-HMM（log-HMMも）のグラフ構造を見直すことにより、以下を実現する。
（a）ガウス混合HMMにおいて、各ガウス（Gaussian）pdfは、到着する状態jに依存する。すべてのN×Kのガウス密度の学習には、多様な長いトレーニング系列を必要とする。
(b)離散型の場合のbjkの役割を連続モードのcjkに割り当てる。また、平均値ベクトルと共分散行列は遷移状態jに依存しないという場合を考慮する。そして、図６の離散型の場合をまっすぐ延長した構造となる。これを半連続HMMと呼んでいる。この構造のもう一つの解釈は、MLE-VQ HMM (Maximum Likelihood Vector Quantization HMM)である。
　半連続のalpha-HMMモデルは、数１４９を変形し、次式となる。
Figure JPOXMLDOC01-appb-M000319
　したがって、半連続のalpha-HMMの更新方程式は以下のようになる。
[単一配列の半連続のalpha-HMMの初期状態確率]
　更新式は、数１３０である。
[単一配列の半連続のalpha-HMMの状態遷移確率]
　更新式は、数１２８である。
[単一配列の半連続のalpha-HMMの分岐確率]
　更新式は、数１５４である。
[単一配列の半連続のalpha-HMMの平均値ベクトル]
　次式のように、数１５６において、μｊｋ｜θl+1（l+1はθの下付きの添え字）をμｊ｜θl+1（l+1はθの下付きの添え字）とし、右辺のkt=kの項を除去したものである。これは、状態ｊに依存しないためである（以後、同様である）。
Figure JPOXMLDOC01-appb-M000320
[単一配列の半連続のalpha-HMMの共分散行列]
　次式のように、数１５８～数１６０において、Σｊｋ｜θl+1（l+1はθの添え字），μｊｋ｜θl（lはθの添え字），μｊｋ｜θl-1（l-1はθの添え字）を、それぞれΣｊ｜θl+1（l+1はθの添え字），μｊ｜θl（lはθの添え字），μｊ｜θl-1（l-1はθの添え字）とし、右辺のkt=kの項を除去したものである。
Figure JPOXMLDOC01-appb-M000321
ここで、各要素は次式である。
Figure JPOXMLDOC01-appb-M000322
Figure JPOXMLDOC01-appb-M000323
　本発明の第５実施例におけるプログラムの処理手順は、前記第３実施例の場合とほぼ同様であるため、図８を用いて説明する。ただし、本実施例では、観測データｙが、連続的な多変数の観測結果として現れる単一配列の半連続系列データである点で異なる。
　本実施例は、前記第３実施例の場合の特例と考えることもできるので、置き換わる数式の対応関係に絞って説明する。まず、繰返し更新される未知パラメータの集合の組み合わせを数１６１とし、推定装置１０が最終的に算出しようとするHMMの確率構造は、数６６の代わりに数１７８とする。ステップＳ３２で決定する未知パラメータの初期値は、数６７の代わりに次式となる。
Figure JPOXMLDOC01-appb-M000324
　ステップＳ３２における収束判定値の決定方法は、数８３と同じである。
　ステップＳ３３では、最初にステップＳ３２で決定した未知パラメータの初期値を利用し、それ以降はステップＳ３９で更新された未知パラメータとカウント値を利用して、次の数１８４に示す各確率値と、数１８５に示す各カウント値を、実際に推定装置１０のメモリ（図示せず）に設定する。
Figure JPOXMLDOC01-appb-M000325
Figure JPOXMLDOC01-appb-M000326
　その後のステップＳ３４～Ｓ３６の演算動作は、前記実施例３と同様に前記実施例１の数７０～数７４の一部の符号を置き換えて適用することによって行なう。
　ステップＳ３７では、l+1回目に繰り返される確率π，a，ｃと、平均値ベクトルμ，及び共分散行列Σの組み合わせを条件として、状態iとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値と、状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂjk（ｙt）において、分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃｊｋ|l+1（l+1はθの添え字）の値と、平均値ベクトルμj|θl+1（l+1はθの添え字）の値と，共分散行列Σj|θl+1（l+1はθの添え字）の値とを、ステップＳ３４で求めたl回目の確率Ｐ（ｙ｜θl）と、それよりも１回前の確率Ｐ（ｙ｜θl-1）を利用し、且つステップＳ３５で求めたl回目のカウント値Ｎａij｜θl，Ｎcjk｜θl（ijはａの添え字、jkはcの添え字、lはθの添え字）と、それよりも１回前のカウント値Ｎａij｜θl-1，Ｎcjk｜θl-1（ijはａの添え字、jkはｃの添え字、l-1はθの添え字）を利用して計算する部分であり、これは前記数１３０、数１２８、数１５４、数１７９～数１８２に対応している。
　このステップＳ３７においても高速化パラメータの値βが用いられ、添え字l+1の部分の確率を計算する上で、時間シフトした添え字lや添え字l-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したこと等は、前記実施例３と同様にその効果として挙げられる。
　本実施例では、配列が１本の観測データをレジスタ１に格納しているが、その場合のステップＳ３７における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数１３０、状態遷移確率の更新式は数１２８、分岐確率の更新式は数１５４、平均値ベクトルの更新式は数１７９、共分散行列の更新式は数１８０～数１８２である。
　ステップＳ３８では、前記実施例３と同様に、前記ステップ３７で計算された新たな確率量に基づく尤度Ｐを用いて、その尤度ＰとステップＳ３２における収束判定値との比較により、ステップＳ３３～Ｓ３９の計算が収束したか否かを判定する。収束してなければ、ステップＳ３９に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように２つの過去値をシフトさせる。
Figure JPOXMLDOC01-appb-M000327
　ここでは、繰り返しの回数が１つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。なお、ステップＳ３３とステップＳ３９を一纏めにして、更新された確率量およびカウント値をそのまま次の繰り返しのために、推定装置１０のメモリに設定記憶させてもよい。
　一方、ステップＳ３８において計算が収束していると判定された場合、ステップＳ４０に移行して、次式に示す計算された５組の確率量を用いたＨＭＭを採用する。
Figure JPOXMLDOC01-appb-M000328
　上記数１８７で計算した各値を用いて、推定装置１０は数１７８に示すＨＭＭの確率構造を出力することができる。
　図８に示すアルゴリズムを実現する推定装置の構成は、図２に示した前記第３実施例の推定装置１０と同様であるため、ここでは説明を省略する。
　以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１４とを備えたＨＭＭの推定装置１０において、推定手段１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを更新設定する更新設定手段２４と、更新設定手段２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数１８４，数１８５を参照）と共に、前記レジスタ１から読み出した観測データと、初期設定手段２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１３０，数１２８，数１５４，数１７９～数１８２を参照）演算手段２６と、演算手段２６による計算の収束を判定し、計算が収束していなければ、演算手段２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段２６で計算した新たな各確率量を最終的な値として出力させる判定手段２８とを備えている。
　このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
　また、ここでの演算手段２６は、l+1回目に状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、前記数１３０の式で計算し、l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（t+1はθの添え字）の値を、l回目およびl-1回目における状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、前記数１２８の式で計算し、状態ｊで状態ｋが出力される出力確率ｂjk｜θl+1（l+1はθの添え字）の際に分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1の値と、平均値ベクトルμjk｜θl+1（l+1はθの添え字）の値と、共分散行列Σjk｜θl+1の値とを、l回目およびl-1回目における前記出力の期待値Ｎcjk｜θlおよびＮcjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、前記数１５４、数１７９～数１８２の式で計算する構成となっている。
　そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
　さらに、ここでの更新設定手段２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
　なお上述した本実施例の作用効果は、初期設定手段２２としての動作を実行する初期設定ステップと、更新設定手段２４としての動作を実行する更新設定ステップと、演算手段２６としての動作を実行する演算ステップと、判定手段２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。
また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。
　本実施例では、半連続的なシンボルで複数配列のデータ系列の場合の実施形態について説明する。
　この場合の複数配列の更新方程式は、平均ベクトルと共分散の状態依存性を制限することによって得られる。初期状態確率と状態遷移確率と分岐確率の更新方程式は、数１４２～数１４５、および数１６７～数１６９と同様である。したがって、次式のようになる。ここで、各記号の上にバー(-)を付し、バー(-）のない場合と記号を区別している。これは、例えば、数１７１と以下の数１８９では数式が異なっており、同じ記号が使えないためである。すなわち、バー(-）それ自体に特別な意味はない。
　また、本実施例では、前記第５実施例において添え字がｊだった記号は、添え字がｋになっている。これは、kt=kの項の代わりにｓｔ＝ｊの項が除去されたためである。（以後、同様である）。
Figure JPOXMLDOC01-appb-M000329
Figure JPOXMLDOC01-appb-M000330
　上記数１８８と同様に、共分散行列の更新方程式は、状態依存性を除去することによって次式のように得られる。
Figure JPOXMLDOC01-appb-M000331
Figure JPOXMLDOC01-appb-M000332
　本発明の第６実施例におけるプログラムの処理手順は、前記第４実施例の場合とほぼ同様であるため、図９を用いて説明する。ただし、本実施例では、観測データｙが、連続的な多変数の観測結果として現れる複数配列の半連続系列データである点で異なる。
　本実施例は、前記第４実施例の場合の特例と考えることもできるので、置き換わる数式の対応関係に絞って説明する。まず、繰返し更新される未知パラメータの集合の組み合わせを数１６１とし、推定装置１００が最終的に算出しようとするHMMの確率構造は、数６６の代わりに数１７８とする。ステップＳ４２で決定する未知パラメータの初期値は、数１８２と同じである。
　ステップＳ４２における収束判定値の決定方法は、数８３と同じである。
　ステップＳ４３では、最初にステップＳ４２で決定した未知パラメータの初期値を利用し、それ以降はステップＳ４９で更新された未知パラメータとカウント値を利用して、次の数１９２に示す各確率値と、数１９３に示す各カウント値を、実際に推定装置１００のメモリ（図示せず）に設定する。
Figure JPOXMLDOC01-appb-M000333
Figure JPOXMLDOC01-appb-M000334
　その後のステップＳ４４～Ｓ４６の演算動作は、前記実施例４と同様に前記実施例２の数８８～数９２の一部の符号を置き換えて適用することによって行なう。
　ステップＳ４７では、l+1回目に繰り返される確率π，a，ｃと、平均値ベクトルμ，及び共分散行列Σの組み合わせを条件として、状態iとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値と、状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値と、状態ｊで状態ｋが出力される出力確率ｂjk（ｙt）において、分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃｊｋ|l+1（l+1はθの添え字）の値と、平均値ベクトルμ-k|θl+1（l+1はθの添え字であり、-はμの上に付く）の値と，共分散行列Σk|θl+1（l+1はθの添え字であり、-はΣの上に付く）の値とを、ステップＳ３４で求めたl回目の確率Ｐ（ｙ（n）｜θl）と、それよりも１回前の確率Ｐ（ｙ（n）｜θl-1）を利用し、且つステップＳ３５で求めたl回目のカウント値Ｎ（n） ａij｜θl，Ｎ（n） cjk｜θl（ijはａの添え字、jkはcの添え字、lはθの添え字）と、それよりも１回前のカウント値Ｎ（n） ａij｜θl-1，Ｎ（n） cjk｜θl-1（ijはａの添え字、jkはｃの添え字、l-1はθの添え字）を利用して計算する部分であり、これは前記数１４２～数１４５、数１６７～数１６９、数１８８～数１９１に対応している。
　このステップＳ４７においても高速化パラメータの値βが用いられ、添え字l+1の部分の確率を計算する上で、時間シフトした添え字lや添え字l-1の値を用いて自己撞着性を取り除くことができたこと、その結果が一つ前の過去値を利用できる計算方法になったこと、その一つ前の過去値の利用が高速性につながったこと、および過去値を利用する重み（高速化パラメータの値）が１≦β＜３に拡張したこと等は、前記実施例４と同様にその効果として挙げられる。
　本実施例では、配列がＭ本の観測データをレジスタ１０１に格納しているが、その場合のステップＳ４７における確率用の更新式は、上述した通りである。すなわち、初期状態確率の更新式は数１４２、状態遷移確率の更新式は数１４３～数１４５、分岐確率の更新式は数１６７～数１６９、平均値ベクトルの更新式は数１８８、数１８９、共分散行列の更新式は数１９０、数１９１である。
　ステップＳ４８では、前記実施例４と同様に、前記ステップＳ４７で計算された新たな確率量に基づく尤度Ｐを用いて、その尤度ＰとステップＳ４２における収束判定値との比較により、ステップＳ４３～Ｓ４９の計算が収束したか否かを判定する。収束してなければ、ステップＳ４９に移行して確率量を含む未知パラメータと事象のカウント値を更新して計算が繰り返される。この際、次式に示すように２つの過去値をシフトさせる。
Figure JPOXMLDOC01-appb-M000335
　ここでは、繰り返しの回数が１つ増えることにより、l回目の確率量およびカウント値がl-1回目の確率量およびカウント値に更新され、l+1回目の確率量およびカウント値がl回目の確率量およびカウント値に更新される。
　一方、ステップＳ４８において計算が収束していると判定された場合、ステップＳ５０に移行して、次式に示す計算された５組の確率量を用いたＨＭＭを採用する。
Figure JPOXMLDOC01-appb-M000336
　上記数１９５で計算した各値を用いて、推定装置１０１は数１７８に示すＨＭＭの確率構造を出力することができる。
　図９に示すアルゴリズムを実現する推定装置の構成は、図５に示した前記第２実施例の推定装置１０１と同様であるため、ここでは説明を省略する。
　以上のように本実施例においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１０１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１１４とを備えたＨＭＭの推定装置１００において、推定手段１１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段１２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ-(-はμの上に付く），共分散行列Σ-(-はΣの上に付く），初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを更新設定する更新設定手段１２４と、更新設定手段１２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いる（数１９２，数１９３を参照）と共に、前記レジスタ１０１から読み出した観測データと、初期設定手段１２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１４２～数１４５，数１６７～数１６９，数１８８～数１９１を参照）演算手段１２６と、演算手段１２６による計算の収束を判定し、計算が収束していなければ、演算手段１２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段１２６で計算した新たな各確率量を最終的な値として出力させる判定手段１２８とを備えている。
　このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ-(-はμの上に付く），共分散行列Σ-(-はμの上に付く），初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
前記数１６７～数１７３の式で計算する構成となっている。
　また、ここでの演算手段１２６は、l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、前記数１４２の式で計算し、l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（t+1はθの添え字）の値を、l回目およびl-1回目における状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、前記数１４３～数１４５の式で計算し、状態ｊで状態ｋが出力される出力確率ｂjk｜θl+1（l+1はθの添え字）の際に分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1の値と、平均値ベクトルμ-k｜θl+1（l+1はθの添え字であり、-はμの上に付く）の値と、共分散行列Σ-k｜θl+1(l+1はθの添え字であり、-はΣの上に付く）の値とを、l回目およびl-1回目における前記分岐の期待値Ｎcjk｜θlおよびＮcjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、前記数１６７～数１６９、数１８８～数１９１の式で計算する構成となっている。
　そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
　さらに、ここでの更新設定手段１２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
　なお上述した本実施例の作用効果は、初期設定手段１２２としての動作を実行する初期設定ステップと、更新設定手段１２４としての動作を実行する更新設定ステップと、演算手段１２６としての動作を実行する演算ステップと、判定手段１２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。
また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。
　上記実施例１～６では、alpha-HMMの６つのタイプ、すなわち、{離散的,連続的,半連続的}×{単一配列,複数配列}について述べた。本実施例は、区切りのある連続アルファベットの場合である。これは、離散的な文字と連続的な文字が混在している場合の解釈である。この離散と連続が混在している更新方程式は、数１５４、数１６９にある総和を２段の総和に変形し、ゼロ出力の許容によって得られる。例えば、数１５４の分母の総和は、次式のような２段の総和に分割される。
Figure JPOXMLDOC01-appb-M000337
　ここで、Dgは、部分が重複しない集合である。Ｇは、その部分の濃度である。もし、Dgが連続的なアルファベットの特定のサブクラスに相当するならば、このことは、そのアルファベットが離散的なシンボルを伴うものと見なすことができる。換言すると、データが連続している各区間において、それぞれ総和を取り、さらに当該総和どうしについて総和を取っている。この離散的で連続的なアルファベットの場合にも、上記実施例のように単一配列の場合と複数配列の場合がある。したがって、この申請は、全部で８つのタイプのalpha-HMMを有している。
　なお、上式の（・）という記号は、Σの中を省略し、総和の範囲を示すΣだけについて考察するために用いたものである。
　この場合には、離散値は連続値を出す集団をグループ化した際のラベルと考えればよい。すなわち、データ系列は次式のように示される。
Figure JPOXMLDOC01-appb-M000338
　ここで、ラベルは、次式として考える。
Figure JPOXMLDOC01-appb-M000339
　上式のＧは、Ｇ≦Ｋ（分岐枝の数）である。そして、新たにｙt-（-はｙtの上に付く）をｙtと書き直せば、ΣT t=1（・）（TはΣの直上に付き、ｔ＝１は直下に付く）がグループ化に対応する部分を含んでいる場合には，これをΣt∈Dg（・）（t∈DgがΣの直下に付く）に変更するだけでよい．ただし、Dg＝｛t|ct＝ｇ｝である。
　以下、各パラメータの更新式について述べる。
・離散連続混在系列の単数配列の場合
[単一配列の離散連続混在系列のalpha-HMMの初期状態確率]
　初期状態確率は数１３０と同じである。
[単一配列の離散連続混在系列のalpha-HMMの状態遷移確率]
　状態遷移確率は数１２８と同じである。
[単一配列の離散連続混在系列のalpha-HMMの分岐確率]
　分岐確率はグループ化情報を併せもち、数１５４に代わって次式となる。以後、これを単にグループ化確率とよぶ。
Figure JPOXMLDOC01-appb-M000340
[単一配列の離散連続混在系列のalpha-HMMの平均値ベクトル]
　各グループに対する平均値ベクトルは数１５６に代わって次式となる。
Figure JPOXMLDOC01-appb-M000341
[単一配列の離散連続混在系列のalpha-HMMの共分散行列]
　各グループに対する共分散行列は、数１５８～数１６０に代わって、次式となる。
Figure JPOXMLDOC01-appb-M000342
　各要素は次式である。
Figure JPOXMLDOC01-appb-M000343
Figure JPOXMLDOC01-appb-M000344
・離散連続混在系列の複数配列の場合
[複数配列の離散連続混在系列のalpha-HMMの初期状態確率]
　初期状態確率は数１４２と同じである。
[複数配列の離散連続混在系列のalpha-HMMの状態遷移確率]
　状態遷移確率は数１４３～数１４５と同じである。
[複数配列の離散連続混在系列のalpha-HMMの分岐確率]
　分岐確率すなわちグループ化確率は、数１６７～数１６９の変形として次式のようになる。
Figure JPOXMLDOC01-appb-M000345
　各要素は次式である。
Figure JPOXMLDOC01-appb-M000346
Figure JPOXMLDOC01-appb-M000347
[複数配列の離散連続混在系列のalpha-HMMの平均値ベクトル]
　平均値ベクトルは、数１７０、数１７１の変形として次式となる。
各要素は次式である。
Figure JPOXMLDOC01-appb-M000349
[複数配列の離散連続混在系列のalpha-HMMの共分散行列]
　共分散行列は、数１７２、数１７３の変形として次式となる。
Figure JPOXMLDOC01-appb-M000350
各要素は次式である。
　本実施例も上記第３～６実施例と同様に、符号や未知パラメータの更新式が異なるだけであり、それらの数式を、単一配列の場合には、図２に示すハードウェア構成と図８に示すフローチャートに適用し、複数配列の場合には、図５に示すハードウェア構成と図９に示すフローチャートに適用するものであるため、ここでは、詳細な説明は省略する。なお、繰返し更新される未知パラメータの集合の組み合わせは数１６１とし、各更新式は、上述したように単一配列の場合は、数１３０、数１２８、数１９９～２０３であり、複数配列の場合は、数１４２～数１４５、数２０４～数２１０である。
　なお、本発明は上記実施例に限定されるものではなく、本発明の趣旨を逸脱しない範囲で変更可能である。例えば、上記各実施例の近似計算において、第３実施例で示した変形例のように遡る過去値を増やしたり、高次の近似にしたりすることも可能である。
　以上のように本実施例（単一配列の場合もほとんど同様であるため、ここでは省略する）においても、入力される観測データを時系列に格納する記憶手段としてのレジスタ１０１と、観測データがどのような確率モデルであるのかを、ＨＭＭの未知パラメータを算出することで推定する推定手段１１４とを備えたＨＭＭの推定装置１００において、推定手段１１４は、ＨＭＭの高速化パラメータの値βを設定する初期設定手段１２２と、前記ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを更新設定する更新設定手段１２４と、更新設定手段１２４で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記レジスタ１０１から読み出した観測データと、初期設定手段１２２で設定した高速化パラメータの値βとを用い、テーラー展開による微小近似を適用して新たな各確率量および各期待値を計算する（数１４２～数１４５，数２０４～数２１０を参照）演算手段１２６と、演算手段１２６による計算の収束を判定し、計算が収束していなければ、演算手段１２６で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、演算手段１２６で計算した新たな各確率量を最終的な値として出力させる判定手段１２８とを備えている。
　このようにすれば、ＨＭＭの未知パラメータとして、状態遷移確率ａ，分岐確率ｃ，平均値ベクトルμ，共分散行列Σ，初期状態確率πおよび尤度Ｐ（ｙ｜θ）の各確率量と、状態遷移の期待値Ｎａおよび分岐の期待値Ｎｃを算出する際に、時間シフトと微小近似を適用して新たな各確率量および各期待値を計算することで、従来の自己撞着の矛盾を回避することができ、特殊な条件下以外であっても、未知パラメータの計算が可能なＨＭＭ推定アルゴリズムを得ることが可能になる。またその形式は、時間シフトした各確率量および各期待値を蓄積された過去情報として利用するだけなので、演算処理の時間を食わず、非常に高速に未知パラメータを求めることができる。
　また、ここでの演算手段１２６は、l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、前記数１４２の式で計算し、l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（t+1はθの添え字）の値を、l回目およびl-1回目における状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、前記数１４３～数１４５の式で計算し、状態jでのグループ化確率ｃjg｜θl+1（l+1はθの添え字）と、そのグループ化を反映した出力確率ｂjg｜θl+1の値と、平均値ベクトルμjg｜θl+1（l+1はθの添え字）の値と、共分散行列Σjg｜θl+1(l+1はθの添え字）の値とを、l回目およびl-1回目における前記分岐の期待値Ｎcjg｜θlおよびＮcjg｜θl-1（jgはｃの添え字、lまたはl-1はθの添え字）を利用して、前記数２０３～数２０９の式で計算する構成となっている。
　そのため、特に高速化パラメータの値βが１でなければならない特殊な条件下以外であっても、未知パラメータの計算が可能になる。
　さらに、ここでの更新設定手段１２４は、高速化パラメータの値を１＜β＜３に設定するのが望ましい。そうすれば、繰り返しの計算が発散しないβが３未満の範囲で、βを１よりも大きく設定して、従来よりも未知パラメータの計算を確実に高速化させることが可能になる。
　なお上述した本実施例の作用効果は、初期設定手段１２２としての動作を実行する初期設定ステップと、更新設定手段１２４としての動作を実行する更新設定ステップと、演算手段１２６としての動作を実行する演算ステップと、判定手段１２８としての動作を実行する判定ステップとを備えたＨＭＭの推定方法であっても、全く同様に発揮されるし、またそうした手段を推定手段１１４として、コンピュータに機能させるＨＭＭの推定プログラムであっても、同様に発揮される。
また、本実施形態は連続系列の観測データにも適用できるため、コンピュータによる音声認識やロボットの動作認識だけでなく、音声の合成やロボットの動作生成等にも応用できる。
　本発明で適用するalpha-ＨＭＭアルゴリズムは、音声認識と合成，ロボットの動作認識と生成，神経情報認識，および生命情報配列の認識などの非常に広範な応用性を有している。具体的には、コンピュータによる音声認識において、観測データとなる音声の特徴パラメータの時間的な変化と確率的な変動とを統計的に扱うために、ここで提案したalpha-ＨＭＭアルゴリズムを使用することができる。
　その中で、alpha-ＨＭＭアルゴリズムを用いた音声認識における学習では、複数の音声を入力としてＨＭＭのパラメータを決定する。そして、実際の音声認識動作のときは、認識対象の音声よりパラメータが決定されたＨＭＭを用いてその確率を計算することが可能になる。また、音声認識等の場合は、離散的なデータに基づいて各確率量を算出するが、連続的なデータに基づいて各確率量を算出すれば、音声の合成や、ロボットの動作の生成が可能となる。
　また、本発明はゲノム配列の認識（例えば、一乃至複数の癌患者からＤＮＡ情報を集めて、それらのデータを上記実施例で提案したalpha-ＨＭＭアルゴリズムに取り込んで癌になる確率を計算する等）、或いはロボットの動作認識などにも適用することができる。

Claims (23)

  1. 　観測データを時系列に格納する記憶手段と、
    　前記観測データがどのような確率モデルであるのかを、隠れマルコフモデルの未知パラメータを算出することで推定する推定手段とを備えた隠れマルコフモデルの推定装置において、
    　前記推定手段は、前記隠れマルコフモデルの高速化パラメータを設定する初期設定手段と、
    　前記隠れマルコフモデルの未知パラメータとして、状態遷移，出力，初期状態および尤度の各確率量と、状態遷移および出力の各期待値を更新設定する更新設定手段と、
    　前記更新設定手段で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記記憶手段から読み出した観測データと、前記初期設定手段で設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算手段と、
    　前記演算手段による計算の収束を判定し、計算が収束していなければ、前記演算手段で計算した新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、前記演算手段で計算した新たな各確率量を最終的な値として出力させる判定手段とを備えた隠れマルコフモデルの推定装置。
  2. 　t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂとの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000001
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
    　前記演算手段は、前記新たな各確率量として、
    　ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000002
     　ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000003
     　状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
    Figure JPOXMLDOC01-appb-M000004
     　ことを特徴とする請求項１記載の隠れマルコフモデルの推定装置。
  3. 　前記更新設定手段は、前記高速化パラメータの値を１＜β＜３に設定することを特徴とする請求項２記載の隠れマルコフモデルの推定装置。
  4. 　記憶手段に観測データを時系列に格納し、
    　前記観測データがどのような確率モデルであるのかを、推定手段が隠れマルコフモデルの未知パラメータを算出することで推定する隠れマルコフモデルの推定方法において、
    　前記未知パラメータの算出は、前記隠れマルコフモデルの高速化パラメータを設定する初期設定ステップと、
    　前記隠れマルコフモデルの未知パラメータとして、状態遷移，出力，初期状態および尤度の各確率量と、状態遷移および出力の各期待値を更新設定する更新設定ステップと、
    　前記更新設定ステップで更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記記憶手段から読み出した観測データと、前記初期設定ステップで設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算ステップと、
    　前記演算ステップによる計算の収束を判定し、計算が収束していなければ、前記演算ステップで計算した新たな各確率量および各期待値を前記更新設定ステップで設定更新させ、計算が収束していれば、前記演算ステップで計算した新たな各確率量を最終的な値として出力させる判定ステップとからなる隠れマルコフモデルの推定方法。
  5. 　t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂとの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000005
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
    　前記演算ステップは、前記新たな各確率量として、
    　ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000006
     　ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000007
     　状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
    Figure JPOXMLDOC01-appb-M000008
     　ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
  6. 　t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂとの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000009
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
    　前記演算ステップは、前記新たな各確率量として、
    　ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000010
     　ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000011
     　状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
    Figure JPOXMLDOC01-appb-M000012
     　ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
  7. 　前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σとの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000013
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
    　前記演算ステップは、前記新たな各確率量として、
    　l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000014
     　l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000015
     　l+1回目の状態ｊで分岐ｋに移行する枝へ遷移する確率を特定する分岐確率cjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎcjk｜θlおよびＮcjk｜θl-1（jkはcの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000016
     　前記分岐確率の確率密度関数の平均値である平均値ベクトルμjk｜θl+1（l+1はθの添え字）を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000017
     　l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σjk|θl+1（l+1はθの添え字）を次の式で計算する
    Figure JPOXMLDOC01-appb-M000018
    Figure JPOXMLDOC01-appb-M000019
    Figure JPOXMLDOC01-appb-M000020
     　ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
  8. 　前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000021
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
    　前記演算ステップは、前記新たな各確率量として、
    　l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000022
     　l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000023
    Figure JPOXMLDOC01-appb-M000024
    Figure JPOXMLDOC01-appb-M000025
     　l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率cjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎcjk｜θlおよびＮcjk｜θl-1（jkはcの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000026
    Figure JPOXMLDOC01-appb-M000027
    Figure JPOXMLDOC01-appb-M000028
     　前記分岐確率の確率密度関数の平均値である平均値ベクトルμjk｜θl+1（l+1はθの添え字）を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000029
    Figure JPOXMLDOC01-appb-M000030
     　l+1回目の状態jで状態jで分岐ｋに移行した時の出力の共分散行列Σjk|θl+1（l+1はθの添え字）を次の式で計算する
    Figure JPOXMLDOC01-appb-M000031
    Figure JPOXMLDOC01-appb-M000032
     　ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
  9. 　前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000033
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
    　前記演算ステップは、前記新たな各確率量として、
    　l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000034
     　l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000035
     　l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率cjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎcjk｜θlおよびＮcjk｜θl-1（jkはcの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000036
     　前記分岐確率の確率密度関数の平均値である平均値ベクトルμj｜θl+1（l+1はθの添え字）を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000037
     　l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σｊ|θl+1（l+1はθの添え字）を次の式で計算する
    Figure JPOXMLDOC01-appb-M000038
    Figure JPOXMLDOC01-appb-M000039
    Figure JPOXMLDOC01-appb-M000040
     　ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
  10. 　前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000041
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
    　前記演算ステップは、前記新たな各確率量として、
    　l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000042
     　l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000043
    Figure JPOXMLDOC01-appb-M000044
    Figure JPOXMLDOC01-appb-M000045
     　l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θlおよびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000046
    Figure JPOXMLDOC01-appb-M000047
    Figure JPOXMLDOC01-appb-M000048
     　前記分岐確率の確率密度関数の平均値である平均値ベクトルμ-k｜θl+1（-はμの直上に付くものであり、l+1はθの添え字）を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000049
    Figure JPOXMLDOC01-appb-M000050
     　l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σ-ｋ|θl+1（-はΣの直上に付くものであり、l+1はθの添え字）を次の式で計算する
    Figure JPOXMLDOC01-appb-M000051
    Figure JPOXMLDOC01-appb-M000052
     　ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
  11. 　前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000053
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
    　前記演算ステップは、前記新たな各確率量として、
    　l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000054
     　l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000055
     　Dgを部分が重複しない集合として、l+1回目の状態jでグループ化gが行なわれるグループ化確率ｃjｇ｜θl+1（l+1はθの添え字）の値を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000056
     　前記分岐確率の確率密度関数の平均値である平均値ベクトルμjｇ｜θl+1（l+1はθの添え字）を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000057
     　l+1回目の状態jでグループ化gが行なわれるときの出力の共分散行列Σｊｇ|θl+1（l+1はθの添え字）を次の式で計算する
    Figure JPOXMLDOC01-appb-M000058
    Figure JPOXMLDOC01-appb-M000059
    Figure JPOXMLDOC01-appb-M000060
     　ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
  12. 　前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000061
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
    　前記演算ステップは、前記新たな各確率量として、
    　l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000062
     　l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000063
    Figure JPOXMLDOC01-appb-M000064
    Figure JPOXMLDOC01-appb-M000065
     　Dgを部分が重複しない集合として、l+1回目の状態jでグループ化gが行なわれる確率の際にg番目のグループへ遷移する確率を特定するグループ化確率ｃjｇ｜θl+1（l+1はθの添え字）の値を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000066
    Figure JPOXMLDOC01-appb-M000067
    Figure JPOXMLDOC01-appb-M000068
     　前記分岐確率の確率密度関数の平均値である平均値ベクトルμjｇ｜θl+1（l+1はθの添え字）を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000069
    Figure JPOXMLDOC01-appb-M000070
     　l+1回目の状態jでグループ化gが行なわれるときの出力の共分散行列Σｊｇ|θl+1（l+1はθの添え字）を次の式で計算する
    Figure JPOXMLDOC01-appb-M000071
    Figure JPOXMLDOC01-appb-M000072
     　ことを特徴とする請求項４記載の隠れマルコフモデルの推定方法。
  13. 　前記更新設定ステップで、前記高速化パラメータの値を１＜β＜３に設定することを特徴とする請求項５～１２記載の隠れマルコフモデルの推定方法。
  14. 　記憶手段に時系列に格納された観測データがどのような確率モデルであるのかを、隠れマルコフモデルの未知パラメータを算出することで推定する推定手段として、コンピュータを機能させる隠れマルコフモデルの推定プログラムにおいて、
    　前記推定手段を、前記隠れマルコフモデルの高速化パラメータを設定する初期設定手段と、
    　前記隠れマルコフモデルの未知パラメータとして、状態遷移，出力，初期状態および尤度の各確率量と、状態遷移および出力の各期待値を更新設定する更新設定手段と、
    　前記更新設定手段で更新設定した直前の各確率量および各期待値のみならず、それより前の時間シフトした各確率量および各期待値を用いると共に、前記記憶手段から読み出した観測データと、前記初期設定手段で設定した高速化パラメータとを用い、微小近似を適用して新たな各確率量および各期待値を計算する演算手段と、
    　前記演算手段による計算の収束を判定し、計算が収束していなければ、前記演算手段による新たな各確率量および各期待値を前記更新設定手段で設定更新させ、計算が収束していれば、前記演算手段による新たな各確率量を最終的な値として確定させる判定手段として機能させることを特徴とする隠れマルコフモデルの推定プログラム。
  15. 　t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂとの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000073
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
    　前記演算手段は、前記新たな各確率量として、
    　ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000074
     　ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000075
     　状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
    Figure JPOXMLDOC01-appb-M000076
     　ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
  16. 　t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、前記出力確率ｂとの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000077
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
    　前記演算手段は、前記新たな各確率量として、
    　ｔ+1回目の状態ｉとなる初期状態確率πｉ｜θt+1（t+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000078
     　ｔ+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記状態遷移の期待値Ｎａij｜θtおよびＮａij｜θt-1（ijはａの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000079
     　状態ｊで状態ｋが出力される出力確率ｂjk｜θt+1（t+1はθの添え字）の値を、t回目およびt-1回目における前記出力の期待値Ｎｂjk｜θtおよびＮｂjk｜θt-1（jkはｂの添え字、tまたはt-1はθの添え字）を利用して、次の式で計算する
    Figure JPOXMLDOC01-appb-M000080
     　ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
  17. 　前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、t回目に繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σとの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000081
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
    　前記演算手段は、前記新たな各確率量として、
    　l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000082
     　l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000083
     　l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θlおよびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000084
     　前記分岐確率の確率密度関数の平均値である平均値ベクトルμjk｜θl+1（l+1はθの添え字）を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000085
     　l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σｊｋ|θl+1（l+1はθの添え字）を次の式で計算する
    Figure JPOXMLDOC01-appb-M000086
    Figure JPOXMLDOC01-appb-M000087
    Figure JPOXMLDOC01-appb-M000088
     　ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
  18. 　前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000089
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
    　前記演算手段は、前記新たな各確率量として、
    　l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000090
     　l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000091
    Figure JPOXMLDOC01-appb-M000092
    Figure JPOXMLDOC01-appb-M000093
     　l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θlおよびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000094
    Figure JPOXMLDOC01-appb-M000095
    Figure JPOXMLDOC01-appb-M000096
     　前記分岐確率の確率密度関数の平均値である平均値ベクトルμjk｜θl+1（l+1はθの添え字）を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000097
    Figure JPOXMLDOC01-appb-M000098
     　l+1回目の状態jで分岐ｋに移行した時の出力の共分散行列Σｊｋ|θl+1（l+1はθの添え字）を次の式で計算する
    Figure JPOXMLDOC01-appb-M000099
    Figure JPOXMLDOC01-appb-M000100
     　ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
  19. 　前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000101
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
    　前記演算手段は、前記新たな各確率量として、
    　l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000102
     　l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000103
     　l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θlおよびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    　前記分岐確率の確率密度関数の平均値である平均値ベクトルμj｜θl+1（l+1はθの添え字）を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000105
     　l+1回目の状態jで分岐ｋに移行した時の出力確率の共分散行列Σｊ|θl+1（l+1はθの添え字）を次の式で計算する
    Figure JPOXMLDOC01-appb-M000106
    Figure JPOXMLDOC01-appb-M000107
    Figure JPOXMLDOC01-appb-M000108
     　ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
  20. 　前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000109
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
    　前記演算手段は、前記新たな各確率量として、
    　l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000110
     　l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000111
    Figure JPOXMLDOC01-appb-M000112
    Figure JPOXMLDOC01-appb-M000113
     　l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率を特定する分岐確率ｃjk｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記分岐の期待値Ｎｃjk｜θlおよびＮｃjk｜θl-1（jkはｃの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000114
    Figure JPOXMLDOC01-appb-M000115
    Figure JPOXMLDOC01-appb-M000116
     　前記分岐確率の確率密度関数の平均値である平均値ベクトルμ-k｜θl+1（-はμの直上に付くものであり、l+1はθの添え字）を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000117
    Figure JPOXMLDOC01-appb-M000118
     　l+1回目の状態ｊで分岐kに移行する枝へ遷移する確率の共分散行列Σ-ｋ|θl+1（-はΣの直上に付くものであり、l+1はθの添え字）を次の式で計算する
    Figure JPOXMLDOC01-appb-M000119
    Figure JPOXMLDOC01-appb-M000120
     　ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
  21. 　前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000121
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとしたときに、
    　前記演算手段は、前記新たな各確率量として、
    　l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000122
     　l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000123
     　Dgを部分が重複しない集合として、l+1回目の状態jでグループ化gが行なわれるグループ化確率ｃjｇ｜θl+1（l+1はθの添え字）の値を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000124
     　前記グループ化確率の確率密度関数の平均値である平均値ベクトルμjｇ｜θl+1（l+1はθの添え字）を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000125
     　l+1回目の状態jでグループ化gが行なわれるときの出力の共分散行列Σｊｇ|θl+1（l+1はθの添え字）を次の式で計算する
    Figure JPOXMLDOC01-appb-M000126
    Figure JPOXMLDOC01-appb-M000127
    Figure JPOXMLDOC01-appb-M000128
     　ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
  22. 　前記未知パラメータとして、さらに分岐確率と平均値ベクトルと共分散行列のパラメータを含めて計算し、繰り返し更新される前記初期状態確率πと、前記状態遷移確率ａと、分岐確率ｃと、平均値ベクトルμと、共分散行列Σの集合の組合せを次の式で示し、
    Figure JPOXMLDOC01-appb-M000129
     　ｓ（太字）を状態ｓの集合とし、ｙ（太字）を前記観測データの値ｙの集合とし、前記高速化パラメータの値をβとし、前記観測データｙ（太字）の配列をＭ本とし、前記Ｍ本の配列のうち、何本目の配列かを示すインデックスをｎとしたときに、
    　前記演算手段は、前記新たな各確率量として、
    　l+1回目の状態ｉとなる初期状態確率πｉ｜θl+1（l+1はθの添え字）の値を、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000130
     　l+1回目の状態ｉから状態ｊに移る状態遷移確率ａij｜θl+1（l+1はθの添え字）の値を、l回目およびl-1回目における前記状態遷移の期待値Ｎａij｜θlおよびＮａij｜θl-1（ijはａの添え字、lまたはl-1はθの添え字）を利用して、次の式で計算し、
    Figure JPOXMLDOC01-appb-M000131
    Figure JPOXMLDOC01-appb-M000132
    Figure JPOXMLDOC01-appb-M000133
     　Dgを部分が重複しない集合として、l+1回目の状態jでグループ化gが行なわれる確率の際にg番目のグループへ遷移する確率を特定するグループ化確率ｃjｇ｜θl+1（l+1はθの添え字）の値を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000134
    Figure JPOXMLDOC01-appb-M000135
    Figure JPOXMLDOC01-appb-M000136
     　前記グループ化確率の確率密度関数の平均値である平均値ベクトルμjｇ｜θl+1（l+1はθの添え字）を次の式で計算し、
    Figure JPOXMLDOC01-appb-M000137
    Figure JPOXMLDOC01-appb-M000138
     　l+1回目の状態jでグループ化gが行なわれるときの出力の共分散行列Σｊｇ|θl+1（l+1はθの添え字）を次の式で計算する
    Figure JPOXMLDOC01-appb-M000139
    Figure JPOXMLDOC01-appb-M000140
     　ことを特徴とする請求項１４記載の隠れマルコフモデルの推定プログラム。
  23. 　前記更新設定手段は、前記高速化パラメータの値を１＜β＜３に設定することを特徴とする請求項１５～２２記載の隠れマルコフモデルの推定プログラム。
PCT/JP2011/058312 2010-07-14 2011-03-31 隠れマルコフモデルの推定方法，推定装置および推定プログラム WO2012008184A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2012524465A JP5709179B2 (ja) 2010-07-14 2011-03-31 隠れマルコフモデルの推定方法，推定装置および推定プログラム

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2010159949 2010-07-14
JP2010-159949 2010-07-14

Publications (1)

Publication Number Publication Date
WO2012008184A1 true WO2012008184A1 (ja) 2012-01-19

Family

ID=45469198

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2011/058312 WO2012008184A1 (ja) 2010-07-14 2011-03-31 隠れマルコフモデルの推定方法，推定装置および推定プログラム

Country Status (2)

Country Link
JP (1) JP5709179B2 (ja)
WO (1) WO2012008184A1 (ja)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111008238A (zh) * 2019-11-15 2020-04-14 武汉楚誉科技股份有限公司 基于关联演化大数据的关键模式自动定位与预警方法
CN111259261A (zh) * 2020-01-02 2020-06-09 中国铁道科学研究院集团有限公司通信信号研究所 基于状态迁移预测的高铁行车网络协同报警优化方法
CN117405573A (zh) * 2023-12-15 2024-01-16 长沙矿冶研究院有限责任公司 一种考虑脉冲数不稳定的射线矿浆浓度计校准方法及系统

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110450154A (zh) * 2019-07-18 2019-11-15 广州弘度信息科技有限公司 一种基于动力学的机械臂柔性控制方法
CN113934930B (zh) * 2021-10-08 2023-04-18 宜兴市旭航电子有限公司 基于概率和隐马尔可夫模型的用户喜好预测系统及方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006285881A (ja) * 2005-04-05 2006-10-19 Sony Corp 学習装置、学習方法、およびプログラム

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3587966B2 (ja) * 1996-09-20 2004-11-10 日本電信電話株式会社 音声認識方法、装置そよびその記憶媒体
US7139703B2 (en) * 2002-04-05 2006-11-21 Microsoft Corporation Method of iterative noise estimation in a recursive framework
JP4705414B2 (ja) * 2005-06-13 2011-06-22 日本電信電話株式会社 音声認識装置、音声認識方法、音声認識プログラムおよび記録媒体
JP2010078650A (ja) * 2008-09-24 2010-04-08 Toshiba Corp 音声認識装置及びその方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006285881A (ja) * 2005-04-05 2006-10-19 Sony Corp 学習装置、学習方法、およびプログラム

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
HIROYUKI SHIOYA ET AL.: "Transformation of Convex Functions and Inequalities of Divergences", THE IEICE TRANSACTIONS ON FUNDAMENTALS OF ELECTRONICS, COMMUNICATIONS AND COMPUTER SCIENCES (JAPANESE EDITION). A, March 1997 (1997-03-01), pages 509 - 515 *
MATSUYAMA Y. ET AL.: "Alpha-EM gives fast Hidden Markov Model estimation: Derivation and evaluation of alpha-HMM, Neural Networks (IJCNN)", THE 2010 INTERNATIONAL JOINT CONFERENCE ON, July 2010 (2010-07-01), pages 1 - 8 *
MATSUYAMA Y. ET AL.: "Fast learning by the a-ECME algorithm, Neural Information Processing", PROCEEDINGS. ICONIP '99. 6TH INTERNATIONAL CONFERENCE, 1999, pages 1184 - 1190 *
MATSUYAMA Y.: "The a-EM algorithm: surrogate likelihood maximization using a-logarithmic information measures", IEEE TRANSACTIONS ON INFORMATION THEORY, March 2003 (2003-03-01), pages 692 - 706 *
MOLINA C.: "Maximum Entropy-Based Reinforcement Learning Using a Confidence Measure in Speech Recognition for Telephone Speech, Audio, Speech, and Language Processing", IEEE TRANSACTIONS ON, July 2010 (2010-07-01), pages 1041 - 1052 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111008238A (zh) * 2019-11-15 2020-04-14 武汉楚誉科技股份有限公司 基于关联演化大数据的关键模式自动定位与预警方法
CN111008238B (zh) * 2019-11-15 2023-10-10 武汉楚誉科技股份有限公司 基于关联演化大数据的关键模式自动定位与预警方法
CN111259261A (zh) * 2020-01-02 2020-06-09 中国铁道科学研究院集团有限公司通信信号研究所 基于状态迁移预测的高铁行车网络协同报警优化方法
CN111259261B (zh) * 2020-01-02 2023-09-26 中国铁道科学研究院集团有限公司通信信号研究所 基于状态迁移预测的高铁行车网络协同报警优化方法
CN117405573A (zh) * 2023-12-15 2024-01-16 长沙矿冶研究院有限责任公司 一种考虑脉冲数不稳定的射线矿浆浓度计校准方法及系统
CN117405573B (zh) * 2023-12-15 2024-03-22 长沙矿冶研究院有限责任公司 一种考虑脉冲数不稳定的射线矿浆浓度计校准方法及系统

Also Published As

Publication number Publication date
JPWO2012008184A1 (ja) 2013-09-05
JP5709179B2 (ja) 2015-04-30

Similar Documents

Publication Publication Date Title
EP3504666B1 (en) Asychronous training of machine learning model
US11797822B2 (en) Neural network having input and hidden layers of equal units
Collins et al. Exponentiated gradient algorithms for conditional random fields and max-margin markov networks
JP4594551B2 (ja) 統合された確率的ランゲージ・モデルを用いたドキュメント画像復号化方法
Papaspiliopoulos et al. Retrospective Markov chain Monte Carlo methods for Dirichlet process hierarchical models
Lavergne et al. Practical very large scale CRFs
Uribe et al. Cross-entropy-based importance sampling with failure-informed dimension reduction for rare event simulation
Cappé Online expectation maximisation
JP5709179B2 (ja) 隠れマルコフモデルの推定方法，推定装置および推定プログラム
Allahverdyan et al. Comparative analysis of viterbi training and maximum likelihood estimation for hmms
Reddy et al. Trust-tech-based expectation maximization for learning finite mixture models
JP2022530447A (ja) ディープラーニングに基づく中国語単語分割方法、装置、記憶媒体及びコンピュータ機器
Cholewa et al. Estimation of the number of states for gesture recognition with Hidden Markov Models based on the number of critical points in time sequence
Jasinska-Kobus et al. Probabilistic label trees for extreme multi-label classification
Sokolovska et al. Efficient learning of sparse conditional random fields for supervised sequence labeling
US20060100874A1 (en) Method for inducing a Hidden Markov Model with a similarity metric
Li Conditional random field
Lifshits et al. Speeding up HMM decoding and training by exploiting sequence repetitions
De Gooijer et al. Kernel-based hidden Markov conditional densities
Winter et al. Machine learning and the future of bayesian computation
Dayar et al. On the numerical analysis of stochastic Lotka-Volterra models
JP2017220001A (ja) 予測装置、予測方法及びプログラム
Mozes et al. Speeding up HMM decoding and training by exploiting sequence repetitions
Tang Autoregressive hidden Markov model with application in an El Nino study
Jaeger et al. Efficient estimation of OOMs

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 11806517

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 2012524465

Country of ref document: JP

122 Ep: pct application non-entry in european phase

Ref document number: 11806517

Country of ref document: EP

Kind code of ref document: A1